CALCULATING ELLIPSE OVERLAP AREAS 



Gary B. Hughes 

California Polytechnic State University 

Statistics Department 
San Luis Obispo, CA 93407-0405, USA 

MOHCINE ChRAIBI 

Jiilich Supercomputing Centre 
Forschungszentrum Jiilich GmbH 
D-52425 Jiilich, Germany 



Abstract. We present a general algorithm for finding the overlap area be- 
tween two ellipses. The algorithm is based on finding a segment area (the 
area between an ellipse and a secant line) given two points on the ellipse. The 
Gauss-Green formula is used to determine the ellipse sector area between two 
points, and a triangular area is added or subtracted to give the segment area. 
For two ellipses, overlap area is calculated by adding the areas of appropriate 
sectors and polygons. Intersection points for two general ellipses are found 
using Ferrari's quartic formula to solve the polynomial that results from com- 
bining the two ellipse equations. All cases for the number of intersection points 
(0, 1, 2, 3, 4) are handled. The algorithm is implemented in c-code, and has 
been tested with a range of input ellipses. The code is efficient enough for use 
in simulations that require many overlap area calculations. 

1. Introduction. Ellipses are useful in many applied scenarios, and in widely dis- 
parate fields. In our research, which happens to be in two very different areas, we 
have encountered a common need for efficiently calculating the overlap area between 
two ellipses. 

In one case, the design for a solar calibrator on-board an orbiting satellite required 
an efficient algorithm for ellipse overlap area. Imaging systems aboard satellites rely 
on semi-conductor detectors whose performance changes over time due to many 
factors. To produce consistent data, some means of calibrating the detectors is 
required; see, e.g., [I]. Some systems use the sun as a light source for calibration. 
In a typical solar calibrator, incident sunlight passes through an attenuator grating 
and impinges on a diffuser plate, which is oriented obliquely to the attenuator 
grating. The attenuator grating is a pattern of circular openings. When sunlight 
passes through the circular openings, projections of the circles onto the oblique 
diffuser plate become small ellipses. The projection of the large circular entrance 
aperture on the oblique diffuser plate is also an ellipse. The total incident light 
on the calibrator is proportional to the sum of all the areas of the smaller ellipses 
that are contained within the larger entrance aperture ellipse. However, as the 
calibration process proceeds, the satellite is moving through its orbit, and the angle 
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from the sun into the calibrator changes (~7° in 2 minutes). The attenuator grating 
elhpses thus move across the entrance aperture, and some of the smaller ellipses 
pass in and out of the entrance aperture ellipse during calibration. Movement of the 
small ellipses across the aperture creates fluctuations in the total amount of incident 
sunlight reaching the calibrator in the range of 0.3 to 0.5%. This jitter creates errors 
in the calibration algorithms. In order to model the jitter, an algorithm is required 
for determining the overlap area of two ellipses. Monte Carlo integration had been 
used; however, the method is numerically intensive because it converges very slowly, 
so it was not an attractive approach for modeling the calibrator due to the large 
number of ellipses that must be modeled. 

In a more down-to-earth setting, populated places such as city streets or build- 
ing corridors can become quite congested while crowds of people are moving about. 
Understanding the dynamics of pedestrian movement in these scenarios can be ben- 
eficial in many ways. Pedestrian dynamics can provide critical input to the design 
of buildings or city infrastructure, for example by predicting the effects of specific 
crowd management strategies, or the behavior of crowds utilizing emergency escape 
routes. Current research in pedestrian dynamics is making steady progress toward 
realistic modeling of local movement; see, e.g., [2]. The model presented in [2] is 
based on the concept of elliptical volume exclusion for individual pedestrians. Each 
model pedestrian is surrounded by an elliptical footprint area that the model uses to 
anticipate obstacles and other pedestrians in or near the intended path. The foot- 
print area is influenced by an individuals' velocity; for example, the exclusion area 
in front of a fast-moving pedestrian is elongated when compared to a slower-moving 
individual, since a pedestrian is generally thinking a few steps ahead. As pedestri- 
ans travel through a confined space, their collective exclusion areas become denser, 
and the areas will eventually begin to overlap. A force-based model will produce a 
repulsive force between overlapping exclusion areas, causing the pedestrians to slow 
down or change course when the exclusion force becomes large. Implementing the 
force-based model with elliptical exclusion areas in a simulation requires calculating 
the overlap area between many different ellipses in the most general orientations. 
The ellipse area overlap algorithm must also be efficient, so as not to bog down the 
simulation. 

Simulations for both the satellite solar calibrator and force-based pedestrian dy- 
namic model require efficient calculation of the overlap area between two ellipses. In 
this paper, we provide an algorithm that has served well for both applications. The 
core component of the overlap area algorithm is based on determining the area of 
an ellipse segment, which is the area between a secant line and the ellipse boundary. 
The segment algorithm forms the basis of an application for calculating the overlap 
area between two general ellipses. 

2. Ellipse area, sector area and segment area. 

2.1. Ellipse Area. Consider an ellipse that is centered at the origin, with its axes 
aligned to the coordinate axes. If the semi-axis length along the a;-axis is A, and 
the semi-axis length along the y-axis is B, then the ellipse is defined by a locus of 
points that satisfy the implicit polynomial equation 
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The same ellipse can be defined parametrically by: 

I 0<t<2. (2) 
y ~ B ■ sm(t) J ^ ' 

The area of such an ellipse can be found using the parameterized form with the 
Gauss-Green formula: 



1 



B 



Area j [x{t) ■ y'{t) - y{t) ■ x'{t)]dt 



1 



A ■ cos(t) • B ■ cos(i) - B ■ sin(t) • {-A) ■ 8m{t)]dt 



(3) 



A-B r^'^ 2, , . 2/ M , A-B f^"" ^ 
cos^(i) + sm^{t)]dt = —— / dt 





--t:-A-B 



2.2. Ellipse Sector Areas. We define the ellipse sector between two points (cci, 
yi) and {x2, 2/2) on the ellipse as the area that is swept out by a vector from the 
origin to the ellipse, beginning at (xi, yi), as the vector travels along the ellipse 
in a counter-clockwise direction from (xi, yi) to (a;2, 2/2)- An example is shown in 
Fig. 1. The Gauss-Green formula can also be used to determine the area of such an 
ellipse sector. 



A ■ B 

Sector Area = — ^ — / dt 



(4) 



3 1 

(.v„>V) 




Figure 1 . The area of an ellipse sector between two points on the 
ellipse is the area swept out by a vector from the origin to the first 
point as the vector travels along the ellipse in a counter-clockwise 
direction to the second point. The area of an ellipse sector can be 
determined with the Gauss-Green formula, using the parametric 
angles 9i and 62- 
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The parametric angle 9 that is formed between the 2:-axis and a point {x, y) on 
the eUipse is found from the eUipse parameterizations: 

x=A-cos{6) =^ 9 = cos^^ {x / A) 

y=B-sm{9) =^ 9 = smT^ {y / B) 

For a circle {A — B \n the ellipse implicit polynomial form) , the parametric angle 
corresponds to the geometric (visual) angle that a line from the origin to the point 
(x, y) makes with the x-axis. However, the same cannot be said for an ellipse; that 
is, the geometric (visual) angle is not the same as the parametric angle used in the 
area calculation. For example, consider the ellipse in Fig. 1; the implicit polynomial 
form is 

5 + fJ = l (5) 
Suppose the point (a;i, j/i) is at (4/a/5, 4/a/5) . The point is on the ellipse, since 

(4/V5)' (4/V5)' _ 475 475 _1 4_ 
42 ^ 22 ~ ^ + ~QP~ "5^5^"^ 

A line segment from the origin to (4/-\/5, 4/a/5) forms an angle with the x-axis 
of 7r/4 (wO. 7485398). However, the ellipse parametric angle to the same point is: 



'4/V5^ 

,V5 



1.10715 



1.10715 



The same angle can also be found from the parametric equation for y: 

„ . _i A/V5\ . -1 / 2 

6^ = sm = sm —p^ 

\ 2 ) VV5y 

The angle found by using the parametric equations does not match the geometric 
angle to the point that defines the angle. 

When determining the parametric angle for a given point (x, y) on the ellipse, 
the angle must be chosen in the proper quadrant, based on the signs of x and y. For 
the ellipse in Fig. 1, suppose the point {x2^ 2/2) is at (—3, — \/7/2). The parametric 
angle that is determined from the equation for x is: 

= cos^i ( — I « 2.41886 



4 

The parametric angle that is determined from the equation for y is: 

.^sin-(z^) ^"""(^) 

The apparent discrepancy is resolved by recalling that inverse trigonometric func- 
tions are usually implemented to return a 'principal value' that is within a conven- 
tional range. The typical (principal-valued) 9 = arccos(a;) function returns angles 
in the range — 9 — 7:, and the typical (principal-valued) 9 — arcsin(a;) function 
returns angles in the range -7r/2 = 9 = 7r/2. When the principal- valued inverse 
trigonometric functions return angles in the typical ranges, the ellipse parametric 
angles, defined to be from the a;-axis, with positive angles in the counter-clockwise 
direction, can be found with the relations in Table 2.2. 
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Quadrant II (a: < and y > 0) 

9 = arccos(a;/A) 
= TT — arcsin(y/i3) 


Quadrant I (a; > and y > 

0) 

9 = arccos(a::/^) 
= arcsin(j//i?) 


Quadrant III (a; < and y < 0) 

9 — 2tt — arccos(a;/A) 
= TT — arcsin(y/i3) 


Quadrant IV (a; > y < 

0) 

9 = 2Tr — arccos(a:/yl) 
= 27r + arcsin(y/i3) 



Table 1. Relations for finding the parametric angle that corre- 
sponds to a given point (a:, y) on the ellipse x'^/A'^ + y^/B^ = 1. 
The parametric angle is formed between the positive a;-axis and a 
line drawn from the origin to the given point, with counterclock- 
wise being positive. For the standard (principal-valued) inverse 
trigonometric functions, the resulting angle will be in the range 
< 9 < 2% for any point on the ellipse. 



The point at (—3, — -\/7/2) on the ellipse of Fig. 1 is in Quadrant III. Using the 
relations in Table 2.2, the parametric angle that is determined from the equation 
for X is: 

9^2-11- arccosf — ) « 3.86433 
4 

The parametric angle that is determined from the equation for y is: 



-77/2 



) w 3.86433 



With the proper angles, the Gauss-Green formula can be used to determine the 
area of the sector from the point at (4/ -\/5, 4/-\/5) to the point (—3, — \/7/2) in the 
ellipse of Fig. 1. 



Sector Area =- 



h)-A-B 



(2n — arccos (^)) 



•4-2 



(6) 



«11.0287 

The Gauss-Green formula is sensitive to the direction of integration. For the 
larger goal of determining ellipse overlap areas, we define the ellipse sector area to 
be calculated from the first point (xi, yi) to the second point {x2, 2/2) in a counter- 
clockwise direction along the ellipse. For example, if the points (a;i, xi) and (a;2, 
7/2) of Fig. 1 were to have their labels switched, then the ellipse sector defined by 
the new points will have an area that is complementary to that of the sector in 
Fig. 1, as shown in Fig. 2. 

Switching the point labels, as shown in Fig. 2, also causes the angle labels to 
be switched, resulting in the condition that 9i > 92- Since using the definitions 
in Table 2.2 will always produce an angle in the range = < 27r for any point 
on the ellipse, the first angle can be transformed by subtracting 2tt to restore the 
condition that 9i < ^2- The sector area formula given above can then be used, 
with the integration angle from {9i - 27t) through 02 . With the angle labels shown 



6 



GARY B. HUGHES AND MOHCINE CHRAIBI 




Figure 2. We define the ellipse sector area to be calculated from 
the first point {xi, yi) to the second point (x2, 2/2) in a counter- 
clockwise direction along the ellipse. 



in Fig. 2, the area of the sector from the point at (— 3,— V7/2) to the point at 
(4/\/5, 4/\/5) in a counter-clockwise direction is: 

(6*2 - {Oi - 27r)) ■ A-B 
Sector Area =^ ^ ^ " 



(27r-arccos(^)) - (arccos(^) -27r ) •4-2 (7) 



«14.1040 

The two sector areas shown in Fig. 1 and Fig. 2 are complementary, in that they 
add to the total ellipse area. Using the angle labels as shown in Fig. 1 for both 
sector areas: 

Total Area.i^^-^^)-^-^ ■ (0i - (02 - 2.)) • A • 



2 

(27r) A-B 

' 2 
=7r • 4 • 2 

^25.1327 



T^-A-B (8) 



2.3. Ellipse Segment Areas. For the overall goal of determining overlap areas 
between ellipses and other curves, a useful measure is the area of what we will call 
an ellipse segment. A secant line drawn between two points on an ellipse partitions 
the ellipse area into two fractions, as shown in Fig. 1 and Fig. 2. We define the 
ellipse segment as the area confined by the secant line and the portion of the ellipse 
from the first point [xi, yi) to the second point (x2, j/2) traversed in a counter- 
clockwise direction. The segment's complement is the second of the two areas that 
are demarcated by the secant line. For the ellipse of Fig. 1, the area of the segment 
defined by the secant hue through the points [xi, yi) and (x2, 1/2) is the area of 
the sector minus the area of the triangle defined by the two points and the ellipse 
center. To find the area of the triangle, suppose that the coordinates for the vertices 
of are known, e.g., as (xi, j/i), {x2, 2/2) and (x^, y^). Then the triangle area can be 
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found by: 



Triangle Area — — 
^ 2 



det 



Xi X2 X3 

yi 2/2 ys 
1 1 1 



(9) 



= 2 • l^^i • (2/2 - 2/3) - X2 ■ (2/1 - 2/3) + X3 ■ (2/1 - 2/2)1 



In the case where one vertex, say (X3, ^3), is at the origin, then the area formula 
for the triangle can be simplified to: 

Triangle Area = ^ ■ \xi ■ y2 ~ X2 ■ yi\ (10) 

For the case depicted in Fig. 1, subtracting the triangle area from the area of the 
ellipse sector area gives the area between the secant line and the ellipse, i.e., the 
area of the ellipse segment counter-clockwise from (xi, yi) to {x2, 2/2): 



Segment Area = 



A- B 1 



2 • Fi • 2/2 



X2-yi\ 



(11) 



For the ellipse of Fig. 1, with the points at (4/V^,4/V5) and (-3,-^7/2), the 
area of the segment defined by the secant line is: 



[2tt — arccos (^)) — arccos { ^^^ ^ 



4-2 



4 -%/7 
V5 ' 



4 

7! 



9.52865 



For the ellipse of Fig. 2, the area of the segment shown is the sector area plus 
the area of the triangle. 

{62- {0i-2tt)) ■ A - B 1 , , 
Segment Area = h - • l^i • 2/2 - 2^2 • 2/i| (12) 

With the points at (-3,-/7/2) and (4/\/5,4/\/5) the area of the segment is: 



{2tt — arccos (^)) — fa 



4/V5 



271 



4-2 



1 

+ 2 



V5 ' 2 



4 

^ --3 



15.60409411 



For the case shown in Fig. 1 and Fig. 2, the sector areas were shown to be 
complementary. The segment areas are also complementary, since the triangle area 
is added to the sector of Fig. 1, but subtracted from the sector of Fig. 2. Using the 
angle labels as shown in Fig. 1 for both sector areas: 



Total Area = 



3i)-A-B 1 
2 



- 7: • l^^i • 2/2 - 2:2 • 2/1I 



-2it))-A-B 1 



=7r-v4-B = 7r-4-2: 



2 

25.1327 



+ 7^-\xi-y2- X2- 2/1I 



(13) 



The key difference between the cases in Fig. 1 and Fig. 2 that requires the area of 
the triangle to be either subtracted from, or added to, the sector area is the size of 
the integration angle. If the integration angle is less than tt, then the triangle area 
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ELLIPSE_SEGMENT Area Algorithm: 

^ _ arccos , j/i > 

^ ^ 1 27r — arccos (xi/^) , yi < 

arccos {x2 /A) , 2/2 > 

2-K — arccos {X2/A) , j/2 < 
1 , Si < 62 

h - 2n, 01 > 02 

3. Area=-'^ +- ^ ^ ' ■ \xi ■ y2 - X2 ■ yi\ 

where: 

the eUipse imphcit polynomial equation is 

^2 -I- B2 — i 

j4 > is the semi-axis length along the i-axis 

i? > is the semi-axis length along the y-axis 

{xi, yi) is the first given point on the ellipse 

{x2, 2/2) is the second given point on the ellipse 

61 and 02 are the parametric angles corresponding to the points 

(xi, yi) and {x2, ^2) 

Table 2. An outline of the ELLIPSE_SEGMENT area algorithm. 



must be subtracted from the sector area to give the segment area. If the integration 
angle is greater than tt, the triangle area must be added to the sector area. 

2.4. A Core Algorithm for Ellipse Segment Area. A generalization of the 
cases given in Fig. 1 and Fig. 2 suggests a robust approach for determining the 
ellipse segment area defined by a secant line drawn between two given points on the 
ellipse. The ellipse is assumed to be centered at the origin, with its axes parallel to 
the coordinate axes. We define the segment area to be demarcated by the secant 
line and the ellipse proceeding counter-clockwise from the first given point (xi, yi) 
to the second given point {x2, 2/2 )• The ELLIPSE_SEGMENT algorithm is outlined 
in Table 2, with pseudo-code presented in List. 1. The ellipse is passed to the 
algorithm by specifying the semi-axes lengths, ^ > and B > 0. The points are 
passed to the algorithm as (xi, xi) and (x2, 2/2), which must be on the ellipse. 

For robustness, the algorithm should avoid divide- by-zero and inverse-trigonometric 
errors, so data checks should be included. The ellipse parameters A and B must be 
greater than zero. A check is provided to determine whether the points are on the 
ellipse, to within some numerical tolerance, s. Since the points can only be checked 
as being on the ellipse to within some numerical tolerance, it may still be possible 
for the x-values to be slightly larger than A, leading to an error when calling the 
inverse trigonometric functions with the argument x/A. In this case, the algorithm 
checks whether the a:- value close to A or -A, that is within a distance that is less 
than the numerical tolerance. If the closeness condition is met, then the algorithm 
assumes that the calling function passed a value that is indeed on the ellipse near 
the point {A, 0) or {-A, 0), so the value of x is nudged back to A or -A to avoid 
any error when calling the inverse trigonometric functions. The core algorithm, 
including all data checks, is shown in List. 1. 
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Listing 1. The ELLIPSE_SEGMENT algorithm is shown for 
calculating the area of a segment defined by the secant line drawn 
between two given points (xi, yi) and (0:2, 2/2) on the ellipse 
X2/A2 +2/2/-B2 = 1. We define the segment area for this algorithm 
to be demarcated by the secant line and the ellipse proceeding 
counter-clockwise from the first given point (xi, yi) to the second 
given point {x2, 2/2 )• 



1 ELLIPSE^EGMENT (A, B, XI, YI , X2 , Y2) 

2 do if (A or B 0) 

3 then return (-1, ERRORJ^LLIPSEJ'ARAMETERS) : DATA CHECK 

4 2222 2222 

5 do if (|X1 /A + YI /B 1| > or |X1 /A + YI /B 1| > ) 

6 then return (-1, ERROR_POINTS3IOT.ON JILLIPSE) : DATA CHECK 

7 do if (I XI I /A > ) 

8 do if |X1| - A > 

9 then return (-1, ERROR JNVERSE.TRIG ) :DATA CHECK 

10 else do if XI < 

11 then XI -A 

12 else XI A 

13 do if (|X2|/A > ) 

14 do if |X2| - A > 

15 then return (-1, ERROR JNVERSE.TRIG ) :DATA CHECK 

16 else do if X2 < 

17 then X2 -A 

18 else X2 A 

19 do if (YI < 0) : ANGLE QUADRANT FORMULA (TABLE 1) 

20 then 1 2 acos (Xl/A) 

21 else 1 aeos (Xl/A) 

22 do if (Y2 < 0) : ANGLE QUADRANT FORMULA (TABLE 1) 

23 then 2 2 acos (X2/A) 

24 else 2 aeos (X2/A) 

25 do if (1 > 2) :MUST START WITH 1 < 2 

26 then 1 1-2 

27 do if ((2 1) > ) : STORE SIGN OF TRIANGLE AREA 

28 then trsgn +1.0 

29 else trsgn +1.0 

30 area 0.5*(A*B*(2 - 1) trsgn*|Xl*Y2 - X2*Y1 | ) 

31 return (area, NORMAL.TERMINATION) 



An implementation of the ELLIPSE_SEGMENT algorithm written in c-code is 
shown in Appendix 4. The code compiles under Cygwin-1.7.7-1, and returns the 
following values for the two test cases presented in Fig. 1 and Fig. 2: 



Listing 2. Return values for the test cases in Fig. 1 and Fig. 2 

32 ec call_es.c ellipse_scgment.c — o call_es . exc 

33 . / e al 1 _e s 

34 Calling e 1 1 i p s e _s e g m e n t . c 

35 Fig. 1: segment area — 9.52864712, return_valuc — 

36 Fig. 2; segment area — 15.60409411, return_value — 

37 sum of ellipse segments — 25.13274123 

38 ellipse area by pi*A*B= 25.13274123 



3. Extending the Core Segment Algorithm to more General Cases. 

3.1. Segment Area for a (Directional) Line through a General Ellipse. 

The core segment algorithm is based on an ellipse that is centered at the origin 
with its axes aligned to the coordinate axes. The algorithm can be extended to 
more general ellipses, such as rotated and/or translated ellipse forms. Start by 
considering the case for a standard ellipse with semi-major axis lengths of A and 
B that is centered at the origin and with its axes aligned with the coordinate 
axes. Suppose that the ellipse is rotated through a counter-clockwise angle (p, and 
that the ellipse is then translated so that its center is at the point (/i, k). The 
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rotated+translated ellipse could then be defined by the set of parameters {A, B, 
h, k, if), with the understanding that the rotation through ip is performed before 
the translation through (/i, k). The approach for extending the core segment area 
algorithm will be to determine analogs on the standard ellipse corresponding to 
any points of intersection between a shape of interest and the general rotated and 
translated ellipse. To identify corresponding points, features of the shape of interest 
are translated by (-/i, -A:), and then rotated by ~Lp. The translated+rotated features 
are used to determine any points of intersection with a similar ellipse that is centered 
at the origin with its axes aligned to the coordinate axes. Then, the core segment 
algorithm can be called with the translated+rotated intersection points. 

Rotation and translation are affine transformations that are also length- and 
area-preserving. In particular, the semi-axis lengths in the general rotated ellipse 
are preserved by both transformations, and corresponding points on the two ellipses 
will demarcate equal partition areas. Fig. 3 illustrates this idea, showing the ellipse 
of Fig. 1 which has been rotated counter-clockwise through an angle if = Stt/S, then 
translated by (/i, k) = (—6, 4-3). 

3 n 




-3 J 



Figure 3. Translation and rotation are affine transformations that 
are also length-and area-preserving. Corresponding points on the 
two ellipses will demarcate equal partition areas. 



< i < 27r 



Suppose that we desire to find the area of the rotated+translated ellipse sector 
defined by the line y = —x, where the line 'direction' travels from lower-right to 
upper- left, as shown in Fig. 3. We describe an approach for finding a segment in a 
rotated+translated ellipse, based on the core ellipse segment algorithm. 

An ellipse that is centered at the origin, with its axes aligned to the coordinate 
axes, is defined parametrically by 

X = A - cos{t) 
y = B ■ sin(i) 

Suppose the ellipse is rotated through an angle ip, with counter-clockwise being 
positive, and that the ellipse is then to be translated to put its center is at the point 
{h, k). Any point (a;, y) on the standard ellipse can be rotated and translated to 
end up in a corresponding location on the new ellipse by using the transformation: 

cos (ip) —sin {ip) 
sin (p) cos (p) 

Rotation and translation of the original standard ellipse does not change the ellipse 
area, or the semi-axis lengths. One important feature of the algorithms presented 
here is that the semi-axis lengths A and B are in the direction of the x- and y-axes, 







Vtr 







X 




h 




y 


+ 


k 



(14) 
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respectively, in the un-rotated (standard) ellipse. In its rotated orientation, the 
semi-axis length A will rarely be oriented horizontally (in fact, for ip = 7r/4, the 
semi-axis length A will be oriented vertically). Regardless of the orientation of the 
rotated+translated ellipse, the algorithms presented here assume that the values of 
A and B passed into the algorithm represent the semi-axis lengths along the x- and 
y-ax.es, respectively, for the corresponding un-rotated, un-translatcd ellipse. The 
angle ip is the amount of counter-clockwise rotation required to put the ellipse into 
its desired location. Specifying a negative value for (fi will rotate the standard ellipse 
through a clockwise angle. The angle ip can be specified in anywhere in the range 
(-8, -|-8); the working angle in the code will be computed from the given angle, 
modulo 27T, to avoid any potential errors (?) when calculating trigonometric values. 
The translation {h, k) is the absolute movement along the coordinate axes of the 
ellipse center to move a standard ellipse into its desired location. Negative values 
of h move the standard ellipse to the left; negative values of k move the standard 
ellipse down. 

To find the area between the given line and the rotated-f translated ellipse, the two 
curve equations can be solved simultaneously to find any points of intersection. But 
instead of searching for the points of intersection with the rotated+translated ellipse, 
it is more efficient to transform the two given points that define the line back through 
the translation {-h, -k) then rotation through -(p. The new line determined by the 
translated+rotated points will pass through the standard ellipse at points that are 
analogous to where the original line intersects the rotated-|-translated ellipse. 

The transformations required to move the given points {xi, yi) and (x2, 2/2) into 
an orientation with respect to a standard ellipse that is analogous to their orientation 
to the given ellipse are the inverse of what it took to rotate-t-translate the ellipse to 
its desired position. The translation is performed first, then the rotation: 

cos(— <^) —sin (—(/?) 
sin(— yi) cos (—(/?) 

Multiplying the vector by the matrix, and simplifying the negative-angle trig func- 
tions gives the following expressions for the translated-|-rotated points: 

Xi„ =cos [ip] ■ [xi -h)-\- sin (tp) • {yi - k) 
Vio =- sin {<p) ■ {xi -h) + cos (cp) • {yt - k) 

The two new points (a;io,yio) and {x2o,y2o) can be used to determine a line, e.g., 
by the point-slope method: 

j/ = j/io + f " I {x-x,,) (16) 

X2o Xiij 

The equation can also be formulated in an alternative way to accommodate cases 
where the translated-|-rotated line is vertical, or nearly so: 

^ = ^io + ^^^^(?y-?Ao) (17) 

2/20 -2/I0 

Points of intersection are found by substituting the line equations into the standard 
ellipse equation, and solving for the remaining variable. For each case, define the 
slope as: 



Xio 











Xi — h 




_ Vi-k _ 



(15) 



X2o - Xlo t/2o - yio 

Then the two substitutions proceed as follows: 
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y =yio + ruy^ ■ {x - xi„) into ^42 + ;b2 = ^ 



(yio +m.y.r ■ 



A2 



52 
2" 



= 1 



(19) 



+ 



(yiof - 2 • mj/cc • xig ■ + {niy^ ■ x-^^f - 



=0 



2 2 

X =xi^ + m^y ■ [y - y^^) into 72+^ = 1 



(a;io 



v{v-yio)) , 

A2 ^ B2 



yl2 52 
1 



^2 + 52 . (^^^) 
S2 



(20) 



+ 

+ 
=0 



2 • {xi^ ■ m^y - {m^yf ■ 2/I0) 
{xiof - 2 • m^y ■ xi„ ■ yig + {m^y ■ yi^f - A^ 



If the translated+rotated line is not vertical, then use the first equation to find 
the a;- values for any points of intersection. If the translated+rotated line is close to 
vertical, then the second equation can be used to find the y-values for any points 
of intersection. Since points of intersection between the; line and the ellipse are 
determined by solving a quadratic equation ax"^ + bx + c, there are three cases to 
consider: 

1. A = 52 _ 4ac < 0: Complex Conjugate Roots (no points of intersection) 

2. A = 52 — 4ac = 0: One Double Real Root (1 point of intersection; line tangent 
to ellipse) 

3. A = 52 _ 4q(c > Q: Two Real Roots (2 points of intersection; line crosses 
ellipse) 

For the first two cases, the segment area will be zero. For the third case, the two 
points of intersection can be sent to the core segment area algorithm. However, to 
enforce a consistency in area measures returned by the core algorithm, the integra- 
tion direction is specified to be from the first point to the second point. As such, 
the ellipse line overlap algorithm should be sensitive to the order that the points 
are passed to the core segment algorithm. We suggest giving the line a 'direction' 
from the first given point on the line to the second. The line 'direction' can then 
be used to determine which is to be the first point of intersection, i.e., the first 
intersection point is where the line enters the ellipse based on what 'direction' the 
line is pointing. The segment area that will be returned from ELLIPSE_SEGMENT 
by passing the line's entry location as the first intersection point is the area within 
the ellipse to the right of the line's path. 
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The approach outhned above for finding the overlap area between a hne and a 
general eUipse is implemented in the ELLIPSE_LINE_OVERLAP algorithm, with 
pseudo-code shown in List. 3. The ellipse is passed to the algorithm by specifying the 
counterclockwise rotation angle ip and the translation (/i, k) that takes a standard 
ellipse and moves it to the desired orientation, along with the semi-axes lengths, 
A > and > 0. The line is passed to the algorithm as two points on the line, (xi, 
yi) and (x2, 2/2)- The 'direction' of the line is taken to be from (xi, t/i) toward (a;2, 
2/2)- Then, the segment area returned from ELLIPSE.SEGMENT will be the area 
within the ellipse to the right of the line's path. 

Listing 3. The ELLIPSEilNE.OVERLAP algorithm is shown 
for calculating the area of a segment in a general ellipse that 
is defined by a given line. The line is considered to have a 
'direction' that runs from the first given point (xi, yi) to the 
second given point {x2, 2/2)- The line 'direction' determines 
the order in which intersection points are passed to the EL- 
LIPSE_SEGMENT algorithm, which will return the area of the 
segment that runs along the ellipse from the first point to the 
second in a counter-clockwise direction. Any routine that calls the 
algorithm ELLIPSEXINE.OVERLAP must be sensitive to the 
order of points that are passed in. 



39 (Area, Code) <- ELLIPSE\.LINE\.OVERLAP ( A, B , H, K, <p ,X1 , YI , X2 , Y2) 



40 


do if (A < or B < 0) 




41 






42 


then return (-1, ERROR^LLIPSE PARAMETERS) :DATA.CIIECK 


43 






44 


do if ( > 2ir) 




45 






46 


then 1^ ^ modulo 27r) : BRING ip INTO 


-27r <^p<2-K (?) 


47 






48 


do if {\Xl\/A> 2tt) 




49 






50 


then XI •(- -A 




51 






52 


XIO ^ cos(ip) * (XI ~ H) + sin(i^) * (YI - K) 




53 






54 


YIO i sin(ip) * (XI - H) + cos(tp) * (YI - K) 




55 






56 


X20 ^ cos(i^) * (X2 - H) + sin(i/p) * (Y2 - K) 




57 






58 


Y20 < sin(i^) * (X2 - H) + cos(ip) * (Y2 - K) 




59 






60 


do if (|X20-X10| > e ) 


:LINE IS NOT VERTICAL 


61 






62 


then m <- (Y20 - Y10)/(X20 - XIO) 


: STORE QUADRATIC COEFFICIENTS 


63 






64 


a ^ (B^ + (A*m)^) / 




65 






66 


b <- (2.0*(Y10*ni m^*X10)) 




67 






68 


c <- (YIO^ - 2.0*m*Y10*X10 + 


(m*X10)^ — B^) 


69 






70 


else if (1 Y20 YIO | > e) 


:LINE IS NOT HORIZONTAL 


71 






72 


then ni <~ (X20 - X10)/(Y20 - 


YIO) : STORE QUADRATIC COEFFS 


73 






74 


a ^ (A^ + (B*m)^)/B^ 




75 






76 


b ^ (2.0*(X10«n m^*Y10)) 


77 






78 


c <-(X10^ - 2.0*m*Y10*X10 + (m*Y10)^ — A^ ) 


79 
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80 




else return (-1, ERROR_LINE_POINTS) : LINE POINTS TOO CLOSE 


81 










82 


discrim <— — 4.0*a*c 




83 










84 


do 


ii' (discrim < 0.0) 




:LINE DOES NOT CROSS ELLIPSE 


85 










86 




then return ( , NO JNTERSECT ) 


87 










88 




else 11' (discrim > 


0.0) 


:TWO INTERSECTION POINTS 


89 










90 




then rootl <— ( — b - 


' s q r t 


(discrim)) / (2.0*a) 


91 










92 




root2 ^ (-b - 


h s q r t 


(discrim)) / (2.0*a) 


93 










94 




else return (0, TANGENT) 


: LINE TANGENT TO ELLIPSE 


95 










96 


do 


if (|X20-X10| > e) 




: ROOTS ARE X-VALUES 


97 










98 




then do if (XIO < 


X20) 


: ORDER PTS SAME AS LINE DIRECTION 


99 










100 




then X 1 -f— 


root 1 




101 










102 




x2 ^ 


r o o t 2 




103 










104 




else xl ■<— 


r o o t 2 




105 










106 




x2 ^ 


rootl 




107 










108 




else do if (YIO < 


Y20) 


: ROOTS ARE Y-VALUES 


109 










110 




then y 1 <— 


root 1 


: ORDER, PTS SAAIE AS LINE DIRECTION 


111 










112 




y2 ^ 


r o o t 2 




113 










114 




else y 1 <— 


r o o t 2 




115 










116 




y2 ^ 


rootl 




117 










118 


(A 


rca,Codc) ^ ELLIPSE^EGMENT (A, B , xl , y 1 , x2 , y2 ) 


119 










120 


do 


if (Code < NORMAL.TERMINATION) 


121 










122 




then return (—1.0, 


Code) 




123 










124 




else return (Area, 


TWO JNTERSECTION.POINTS ) 



An implementation of the ELLIPSEilNE.OVERLAP algorithm in c-code is 
shown in Appendix 5. The code compiles under Cygwin-1.7.7-1, and returns the 
following values for the test cases presented above in Fig. 3, with both line 'direc- 
tions': 

Listing 4. Return values for the test cases in Fig. 3. 



125 cc call_el.e ellipse_line_overlap.c ellipse_segment.c — o call_el.exe 
126 

127 . / c all _e 1 
128 

129 Calling c 1 1 i p s e _ 1 i n e _ o v c r 1 a p . c 
130 

131 area = 4.07186819, roturn.value = 102 

132 

133 reverse: area = 21.06087304, roturn.value = 102 

134 

135 sum of ellipse segments — 25.13274123 
136 

137 total ellipse area by pi*A*B = 25.13274123 



3.2. Ellipse-Ellipse Overlap Area. The method described above for determining 
the area between a line and an ellipse can be extended to the task of finding the 
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overlap area between two general ellipses. Suppose the two ellipses are defined 
by their semi-axis lengths, center locations and axis rotation angles. Let the two 
sets of parameters {Ai, Bi, hi, ki, ifi) and {A2, B2, /i2, ^2, V2) define the two 
ellipses for which overlap area is sought. The approach presented here will be to 
first translate both ellipses by an amount {-hi , -ki ) that puts the center of the first 
ellipse at the origin. Then, both translated ellipses are rotated about the origin 
by an angle -(pi that aligns the axes of the first ellipse with the coordinate axes; 
see Fig. 4. Intersection points are found for the two translated+rotated ellipses, 
using Ferrari's quartic formula. Finally, the segment algorithm described above is 
employed to find all the pieces of the overlap area. 





8 
7 

>. 6 

/ \ . 5 ■ 

v \ 3 








1 

(. 1 . , . lO. 






10 -9 - 


-7 -6 -5 -4\3 -2 -1 ( 


1 2 


/4 5 




:;j 







Figure 4. Intersection points on each curve are used with the 
ellipse segment area algorithm to determine overlap area, by cal- 
culating the area of appropriate segments, and polygons in certain 
cases. For the case of two intersection points, as shown above, the 
overlap area can be found by adding two segments, as shown in 
Fig. 5. 

For example, consider a case of two general ellipses with two (non-tangential) 
points of intersection, as shown in Fig. 4. The translation-|-rotation transformations 
that put the first ellipse at the origin and aligned with the coordinate axes do not 
alter the overlap area. In the case shown in Fig. 4, the overlap area consists of one 
segment from the first ellipse and one segment from the second ellipse. The segment 
algorithm presented above can be used directly for ellipses centered at the origin 
and aligned with the coordinate axes. As such, the desired segment from the first 
ellipse can be found immediately with the segment algorithm, based on the points 
of intersection. To find the desired segment of the second ellipse, the approach 
presented here further translates and rotates the second ellipse so that the segment 
algorithm can also be used directly. The overlap area for the case shown in Fig. 4 
is equal to the sum of the two segment areas, as shown in Fig. 5. Other cases, 
e.g. with 3 and 4 points of intersection, can also be handled using the segment 
algorithm. 

The overlap area algorithm presented here finds the area of appropriate sector(s) 
of each ellipse, which are demarcated by any points of intersection between the 
two ellipse curves. To find intersection points, the two ellipse equations are solved 
simultaneously. This step can be accomplished by using the implicit polynomial 
forms for each ellipse. The first ellipse equation, in its translated-f-rotated position 
is written as an implicit polynomial using the appropriate semi-axis lengths: 
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Figure 5. The area of overlap between two intersecting ellipses 
can be found by using the ellipse sector algorithm. In the case of 
two (non-tangential) intersection points, the overlap area is equal 
to the sum of two ellipse sectors. The sector in each ellipse is 
demarcated by the intersection points. 



In a general form of this problem, the translation+rotation that puts the first 
ellipse centered at the origin and oriented with the coordinate axes will typically 
leave the second ellipse displaced and rotated. The implicit polynomial form for 
a more general ellipse that is rotated and/or translated away from the origin is 
written in the conventional way as: 

AA ■ + BB ■ X ■ y + CC ■ + DD ■ X + EE ■ y + FF = (22) 

Any points of intersection for the two ellipses will satisfy these two equations 
simultaneously. An intermediate goal is to find the implicit polynomial coefficients 
in Ellipse Eq. 22 that describe the second ellipse after the translation+rotation that 
puts the first ellipse centered at the origin and oriented with the coordinate axes. 
The parameters that describe the second ellipse after the translation+rotation can 
be determined from the original parameters for the two ellipses. The first step is 
to translate the second ellipse center (/12, through an amount {-hi, -fci), then 
rotate the center-point through -(pi to give a new center point (/i2Tfl, k^TR)'- 

h2TR =cos {-tpi) ■ (/i2 - hi) - sin {-ipi) ■ {k2 - ki) 

k2TR=sm{-ipi) ■ {h2- hi) +cos(-</7i) -(^2-^1) 

The coordinates for a point {xtr, yxR) on the second ellipse in its new trans- 
lated+rotated position can be found from the following parametric equations, based 
on an ellipse with semi-axis lengths A2 and B2 that is centered at the origin, then 
rotated and translated to the desired position: 



< i < 27r 



xtr = A2 ■ COS (t) ■ COS {ip2 - <(5i) - -62 • sin (t) • sin ((^2 - Vi) + ^2ts 
yTR = A2 ■ cos [t) ■ sin ((/?2 - <^i) + B2 ■ sin (<) • cos {ip2 - ipi) + A:2Tn 

To find the implicit polynomial coefficients from the parametric form, further 
transform the locus of points {xtr, Vtr) so that they lie on the ellipse {A2, B2, 0, 
0, 0), which is accomplished by first translating (xtr, Vtr) through [-{hi - /12), 
-{ki - k2)) and then rotating the point through the angle -{^pi - ^2)' 
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X =cos {ip2 - (pi) ■ {xTR^ - {hi - h2)) - sin {^p2 - 'Pi) ■ [vtr - {ki - fc2)) 
y =sin {(p2 ~ ipi) ■ {xTR - {hi ~ ft.2)) + cos {ip2 - Pi) ■ {vtr - {h - k2)) 

The locus of points (a;, y) should satisfy the standard ellipse equation with the 
appropriate semi-axis lengths: 

Finally, the implicit polynomial coefficients for Ellipse Eq. 22 are found by sub- 
stituting the expressions for the point (x, y) into the standard ellipse equation, 
yielding the following ellipse equation: 



[cos (v?2 - 'Pi) ■ {xTR - {hi - h2)) - sin {ip2 - Pi) ■ {vtr - {h - k2)) ] 

Al 

[sin {p2 - pi) ■ {xTR - {hi - /i-2)) + cos {p2 - pi) ■ {vtr - {ki - k2))f (24) 

Bl 

=1 

where {xxtr, yxR) are defined as above. Expanding the terms, and then re- 
arranging the order to isolate like terms yields the following expressions for the 
implicit polynomial coefficients of a general ellipse with the set of parameters {A2 , 

B2, h2TR, k2TR, P2 ~ Pi)- 



AA 
BB 
CC 
DD 

EE 



COS^ ((^2 - (yfl) , Sin^ ((^2 - </5l) 



Al Bl 



2 • sin {p2 - pi) ■ cos {p2 - pi) 2 • sin {p2 - pi) ■ cos {p2 - pi) 



Al Bl 



SlT?{Lp2-pi) , C0i?{p2-pi) 



Al Bl 



—2 • cos {p2 


- Pi) ■ \h2TB ■ cos (((92 


- Pi) ^ 


■f fc2TR 


■ sin {p2 


~Pi) ] 




Al 










2 • sa\{p2 - 


Pi) ■ [k2TR ■ cos ((^2 - 


Pi) - 




sin {p2 - 






Bl 










-2 • sin {ip2 


- Pi) ■ [h2TR ■ COS (.^2 


- Pi) - 


f k2TR 


■ sin{p2 


-Pi) ] 




Al 










2 • cos {p2 — 


pi) ■ [h2TR ■ sm{p2 - 


pi) - 


k2TR 


COS {p2 - 


-fi)] 



(25) 



Bl 

_ [fe2rfl • COS {P2 -Pl) + fc2Ta ' {^2 - Pi) f 
^2 

[ft.2Tj;. ■ sin {P2 - Pi) - fc2rfl ' COS {p2 - Pi) f _ . 

For the area overlap algorithm presented in this paper, the points of intersection 
between the two general ellipses are found by solving simultaneously the two implicit 
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polynomials denoted above as Ellipse Eq. 21 and Ellipse Eq. 22. Solving for x in 
the first equation: 

J + ^ = l =^ --±W^Ml-^) (26) 

Substituting these expressions for x into Ellipse Eq. 22 and then collecting terms 
yields a quartic polynomial in y. It turns out that substituting either the positive 
or the negative root gives the same quartic polynomial coefficients, which are: 




where: 

cy [4] ^ 

cy [3] 



cy [4] -y^ + cy [3] ■ y' + cy [2] .y^ + cy[l]-y + cy [0] = (27) 



■ AA^ + B'l ■ [Al ■ {BB^ -2-AA- CC) + B^ ■ CC^] 



cy [2 



=2-Bi- [Bl ■ CC ■ EE + Aj ■ {BB -DD-AA- EE)] 



— =Al ■ \ [Bl ■ (2 ■ AA ■ CC - BB^) + DD^ - 2 ■ AA ■ FF] - 2 ■ Af ■ AaH 
Bl L J 

+BI ■{2-CC-FF + EE^) 

=2 • Si • [Al ■ {AA ■ EE- BB- DD) + EE ■ FF] 

Bl 

cy [0 

Bl 

(28) 



= [Al ■ {Al ■ AA - DD) + FF] ■ [Ai ■ {Ai ■ AA + DD) + FF] 



In theory, the quartic polynomial will have real roots if and only if the two curves 
intersect. If the ellipses do not intersect, then the quartic will have only complex 
roots. Furthermore, any real roots of the quartic polynomial will represent values 
of intersection points between the two ellipse curves. As with the quadratic equation 
that arises in the ellipse- line overlap calculation, the ellipse-ellipse overlap algorithm 
should handle all possible cases for the types of quartic polynomial roots: 

1. Four real roots (distinct or not); the ellipse curves intersect. 

2. Two real roots (distinct or not) and one complex-conjugate pair; the ellipse 
curves intersect. 

3. No real roots (two complex-conjugate pairs); the ellipse curves do not inter- 
sect. 

For the method we present here, polynomial roots are found using Ferrari's quar- 
tic formula. A numerical implementation of Ferrari's formula is given in [3]. For 
complex roots are returned, and any roots whose imaginary part is returned as zero 
is a real root. 

When the polynomial coefficients are constructed as shown above, the general 
case of two distinct ellipses typically results in a quartic polynomial, i.e., the coeffi- 
cient cy[4] is non-zero. However, certain cases lead to polynomials of lesser degree. 
Fortunately, the solver in [3] is conveniently modular, providing separate functions 
BIQUADROOTS, CUBICROOTS and QUADROOTS to handle all the possible 
polynomial cases that arise when seeking points of intersection for two ellipses. 
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If the polynomial solver returns no real roots to the polynomial, then the ellipse 
curves do not intersect. It follows that the two ellipse areas are either disjoint, or 
one ellipse area is fully contained inside the other; all three possibilities are shown 
in Fig. 6. Each sub-case in Fig. 6 requires a different overlap-area calculation, i.e. 
either the overlap area is zero (Case 0-3), or the overlap is the area of the first 
ellipse (Case 0-2), or the overlap is the area of the second ellipse (Case 0-1). When 
the polynomial has no real roots, geometry can be used to determine which specific 
sub-case of Fig. 6 is represented. An efhcient logic starts by determining the relative 
size of the two ellipses, e.g., by comparing the product of semi-axis lengths for each 
ellipse. The area of an ellipse is proportional to the product of its two semi-axis 
lengths, so the relative size of two ellipses can be determined by comparing the 
product of semi-axis lengths: 

{TT-Ai-Bi)a{TT-A2-B2) =^ {Ai-Bi)a{A2-B2), a e {'<',' >'} (29) 

Suppose the first ellipse is larger than the second ellipse, then AiBi > ^2^2- In 
this case, if the second ellipse center (/i2Tfl, ^2Tfl) is inside the first ellipse, then 
the second ellipse is wholly contained within the first ellipse (Case 0-1); otherwise, 
the ellipses are disjoint (Case 0-3). The logic relies on the fact that there are 
no intersection points, which is indicated whenever there are no real solutions to 
the quartic polynomial. To test whether the second ellipse center (/i2Tfi, fer/j) is 
inside the first ellipse, evaluate the first ellipse equation at the point x = /i2Tfl, and 
y = k2TR', if the result is less than one, then the point (/i2Tfi, fc2Tfl) is inside the 
first ellipse. The complete logic for determining overlap area when AiBi > A2B2 
is: 

If the polynomial has no real roots, and AiBi > A2B2, and + ^ir" < 1; then 
the first ellipse wholly contains the second, otherwise the two ellipses are disjoint. 



-1 1 3 




Figure 6. When the quartic polynomial has no real roots, the 
ellipse curves do not intersect. It follows that either one ellipse is 
fully contained within the other, or the ellipse areas are completely 
disjoint, resulting in three distinct cases for overlap area. 

Alternatively, suppose that the second ellipse is larger than the first ellipse, then 
AiBi < ^2^2. If the first ellipse center (0, 0) is inside the second ellipse, then the 
first ellipse is wholly contained within the second ellipse (Case 0-2); otherwise the 
ellipses are disjoint (Case 0-3). Again, the logic relies on the fact that there are no 
intersection points. To test whether (0, 0) is inside the second ellipse, evaluate the 
second ellipse equation at the origin; if the result is less than zero, then the origin 
is inside the second ellipse. The complete logic for determining overlap area when 
AiBi < A2B2 is: 
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If the polynomial has no real roots, and AiBi < A2B2, and FF < 0, then the 
second ellipse wholly contains the first, otherwise the two ellipses are disjoint. 

Suppose that the two ellipses are the same size, i.e., AiBi = In this 

case, when no intersection points exist, the ellipses must be disjoint (Case 0-3). It 
also turns out that the polynomial solver of [3] will return no real solutions if the 
ellipses are identical. This special case is also handled in the overlap area algorithm 
presented below. Pseudo-code for a function NOINTPTS that determines overlap 
area for the cases depicted in Fig. G is shown in Fig. 14. 

If the polynomial solver returns either two or four real roots to the quartic equa- 
tion, then the ellipse curves intersect. For the algorithm presented here, all of the 
various possibilities for the number and type of real roots are addressed by creating 
a list of distinct real roots. The first step is to loop through the entire array of 
complex roots returned by the polynomial solver, and retrieve only real roots, i.e., 
only those roots whose imaginary component is zero. The algorithm presented here 
then sorts the real roots, allowing for an efficient check for multiple roots. As the 
sorted list of real roots is traversed, any root that is 'identical' to the previous root 
can be skipped. 

Each distinct real root of the polynomial represents a y-value where the two 
ellipses intersect. Each y-value can represent either one or two potential points 
of intersection. In the first case, suppose that the polynomial root is y = Si (or 
y = —Bi), then the y-value produces a single intersection point, which is at (0, Bi) 
(or (0, -Bi j). In the second case, if the y-value is in the open interval (—Si, Bi), 
then there are two potential intersection points where the j/-value is on the first 
ellipse: 



Each potential intersection point (xj, yj) is evaluated in the second ellipse equa- 
tion: 

AA-x'f+BB- X., ■ Vi + CC ■ yl + DD ■ + EE ■ y, + FF, i = 1,2 

If the expression evaluates to zero, then the point (x, y) is on both ellipses, i.e., 
it is an intersection point. By checking all points (x, y) for each value of y that 
is a root of the polynomial, a list of distinct intersection points is generated. The 
number of distinct intersection points must be either 0, 1, 2, 3 or 4. The case of 
zero intersection points is described above, with all possible sub-cases illustrated in 
Fig. 6. If there is only one distinct intersection point, then the two ellipses must be 
tangent at that point. The three possibilities for a single tangent point are shown 
in Fig. 7. 

For the purpose of determining overlap area, the cases of or 1 intersection points 
can be handled in the same way. When two intersection points exist, there are three 
possible sub-cases, shown in Fig. 8. It is possible that both of the intersection points 
are tangents (Case 2-1 and Case 2-2). In both of these sub-cases, one ellipse must 
be fully contained within the other. The only other possibility for two intersection 
points is a partial overlap (Case 2-3). 
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Figure 7. When only one intersection point exists, the ehipses 
must be tangent at the intersection point. As with the case of zero 
intersection points, either one elhpse is fully contained within the 
other, or the ellipse areas are disjoint. The algorithm for finding 
overlap area in the case of zero intersection points can also be used 
when there is a single intersection point. 



Figure 8. When two intersection points exist, either both of the 
points are tangents, or the ellipse curves cross at both points. For 
two tangent points, one ellipse must be fully contained within the 
other. For two crossing points, a partial overlap must exist 

Each sub-case in Fig. 8 requires a different overlap-area calculation. When two 
intersection points exist, either both of the points are tangents, or the ellipse curves 
cross at both points. Specifically, when there are two intersection points, if one point 
is a tangent, then both points must be tangents. And, if one point is not a tangent, 
then neither point is a tangent. So, it suffices to check one of the intersection points 
for tangency. Suppose the ellipses are tangent at an intersection point; then, points 
that lie along the first ellipse on either side of the intersection will lie in the same 
region of the second ellipse (inside or outside). That is, if two points are chosen 
that lie on the first ellipse, one on each side of the intersection, then both points will 
either be inside the second ellipse, or they will both be outside the second ellipse. 
If the ellipse curves cross at the intersection point, then the two chosen points will 
be in different regions of the second ellipse. 

A logic based on testing points that are adjacent to a tangent point can be 
implemented numerically to test whether an intersection point is a tangent or a 
cross-point. Starting with an intersection point (x, y), calculate the parametric 
angle on the first ellipse, by the rules in Table 2.2: 
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larccos(a;/Ai) y>0 ^^^^ 

I 27r — arccos(a;/Ai) y<0 

A small perturbation angle is then calculated. For the method presented here, 
we seek to establish an angle that corresponds to a point on the first ellipse that is 
a given distance, approximately 2EPS, away from the intersection point: 

(2 • EPS \ 
, , (31) 

The angle EPS Radian is then used with the parametric form of the first ellipse to 
determine two points adjacent to (x, y): 



Xl ^Ai ■ COs{e + i;PS'Radia„) 

yi =Bi • Sin(6l + -BPS'Radian) , , 

(32) 

X2 • C0S(6I - -BPS'Radian) 

1/2 =Bi ■ Sin{0 - EPSYUdie^n) 

Each of the points is then evaluated in the second ellipse equation: 
test, = AA ■ x\ +BB-x^- y, + CC ■ yi + DD ■ x, + EE ■ y, + FF, i^l,2 (33) 

If the value of test.; is positive, then the point (xi, yi) is outside the second 
ellipse. It follows that the product of the two test-point evaluations testitest2 will 
be positive if the intersection point is a tangent, since at a tangent point both 
test points will be on the same side of the ellipse. The product of the test-point 
evaluations will be negative if the two ellipse curves cross at the intersection point, 
since the test points will be on opposite sides of the ellipse. The function ISTANPT 
implements this logic to check whether an intersection point is a tangent or a cross- 
point; pseudo-code is shown in Fig. 18. 

When there are two intersection points, the ISTANPT function can be used to 
differentiate the case 2-3 (Fig. 8) from the cases 2-1 and 2-2. Either of the two 
known intersection points can be checked with ISTANPT. If the intersection point 
is a tangent, then both of the intersection points must be tangents, so the case is 
either 2-1 or 2-2, and one ellipse must be fully contained within the other. For cases 
2-1 and 2-2, the geometric logic used for or 1 intersection points can also be used, 
i.e., the function NOINTPTS can be used to determine the overlap area for these 
cases. If the two ellipse curves cross at the tested intersection point, then the case 
must be 2-3, representing a partial overlap between the two ellipse areas. 

For case 2-3, with partial overlap between the two ellipses, the approach for find- 
ing overlap area is based on using the two points {xi, yi) and {x2, 2/2) with segment 
the algorithm (Table 2; Fig. 2) to determine the partial overlap area contributed 
by each ellipse. The total overlap area is the sum of the two segment areas. The 
two intersection points divide each ellipse into two segment areas (see Fig. 5). Only 
one sector area from each ellipse contributes to the overlap area. The segment al- 
gorithm returns the area between the secant line and the portion of the ellipse from 
the first point to the second point traversed in a counter-clockwise direction. For 
the overlap area calculation, the two points must be passed to the segment algo- 
rithm in the order that will return the correct segment area. The default order is 
counter-clockwise from the first point (xi, yi) to the second point {x2, j/2)- A check 
is made to determine whether this order will return the desired segment area. First, 
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the parametric angles corresponding to (xi, yi) and {x2, 1/2) on the first ellipse are 
determined, by the rules in Table 2.2: 

6I1 = |arccos(a;i/Ai) yi > ^^^^ 

1 27r — arccos(a;i/v4i) yi < 



72 



|arccos(a;2/Ai) 2/2 > ^^^^ 

1 27r — arccos(a;2/Ai) 2/2 < 

Then, a point between (xi, yi) and (0:2, ^2) that is on the first ellipse is found: 



71 P2 



^ ' (36) 

'1 ^- P2 



a^mid =^1 ■ COS 

2/mid =-Bi • sin 

The point (xmid, j/mid) is on the first ellipse between (xi, yi) and (a;2, 2/2) when 
travelling counter- clockwise from {xi, j/i) and {x2, 2/2)- If (a;mid, 2/mid) is inside 
the second ellipse, then the desired segment of the first ellipse contains the point 
{xmid, 2/mid)- In this casc, the segment algorithm should integrate in the default 
order, counterclockwise from (xi,yi) to {x2, 2/2)- Otherwise, the order of the points 
should be reversed before calling the segment algorithm, causing it to integrate 
counterclockwise from {x2, 2/2) to {xi, yi). The area returned by the segment 
algorithm is the area contributed by the first ellipse to the partial overlap. 

The desired segment from the second ellipse is found in a manner to the first 
ellipse segment. A slight difference in the approach is required because the segment 
algorithm is implemented for ellipses that are centered at the origin and oriented 
with the coordinate axes; but, in the general case the intersection points (xi, yi) 
and {x2, 2/2) lie on the second ellipse that is in a displaced and rotated location. 
The approach presented here translates and rotates the second ellipse to the origin 
so that the segment algorithm can be used. It suffices to translate then rotate the 
two intersection points by amounts that put the second ellipse centered at the origin 
and oriented with the coordinate axes: 

a;iTR =i^i - h2TR) ■ cos((/3i - 1P2) + [yi - A:2Tr) • sin((^2 - fi) 

2/iTR ={xi - /i2Tr) • sin((^i - (^2) + (2/1 - fc2TR) • cos((^i - (P2) , ^ 

(37) 

a;2TR ^{X2 - /i2TR) ' COs((^i - (^2) + (2/2 - k2TVt) ' sin((^2 - -i^l) 
2/2TR =(a;2 - /i2Tr) ■ Sm{ipi - ip2) + (2/2 - fc2TR) • C0S(V51 - ip2) 

The new points (a;iTR, 2/itr) and (a::2TR, 2/2Tr) lie on the second ellipse after a 
translation+rotation that puts the second ellipse at the origin, oriented with the 
coordinate axes. The new points can be used as inputs to the segment algorithm 
to determine the overlap area contributed by the second ellipse. As with the first 
ellipse, the order of the points must be determined so that the segment algorithm 
returns the appropriate area. The default order is counter-clockwise from the first 
point (xiTR, 2/itr) to the second point (a;2TR, 2/2Tr)- A check is made to determine 
whether this order will return the desired segment area. First, the parametric angles 
corresponding to points (xitr, 2/itr) and {x2tr, 2/2Tr) on the second ellipse are 
determined, by the rules in Table 2.2: 
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6- 



6 




arccos(a;iTRM2) ?/itr > 

27r - arccos(a;iTRM2) 2/itr < 



arccos(a:2TRM2) 



2/2TR > 



(38) 



(39) 



[27r - arccos(x2TRM2) 2/2TR < 

Then, a point on the second ellipse between (xitr, ?/itr) and (x2tRi y2TR) is 
found: 



The point (xmid, J/mid) is on the second ellipse between (xitr, 2/itr) and (a;2TR, 
2/2Tr) when travelling counter- clockwise from (xitr, ?/itr) and (a;2TR, y2TR)- The 
new point [xmid, ymid) lies on the centered second ellipse. To determine the desired 
segment of the second ellipse, the new point (Xmid, ?/mid) must be rotated then 
translated back to a corresponding position on the once-translated+rotated second 
ellipse: 



If (a^midfiTj ymidRx) is iusidc thc first ellipse, then the desired segment of the 
second ellipse contains the point (xmid, 2/mid)- In this case, the segment algorithm 
should integrate in the default order, counterclockwise from (a;iTR, J/itr) to (x2tr, 
J/2Tr)- Otherwise, the order of the points should be reversed before calling the 
segment algorithm, causing it to integrate counterclockwise from (a:2TRj y2TR) to 
(a^iTR, J/itr)- The area returned by the segment algorithm is the area contributed 
by the second ellipse to the partial overlap. The sum of the segment areas from the 
two ellipses is then equal to the ellipse overlap area. The TWOINTPTS function 
calculates the overlap area for partial overlap with two intersection points (Case 
2-3); pseudo-code is shown in Fig. 15. 

There are two possible sub-cases for three intersection points, shown in Fig. 9. 
One of the three points must be a tangent point, and the ellipses must cross at the 
other two points. The cases are distinct only in the sense that the tangent point 
occurs with ellipse 2 on the interior side of ellipse 1 (Case 3-1), or with ellipse 2 on 
the exterior side of ellipse 1 (Case 3-2). The overlap area calculation is performed 
in the same manner for both cases, by calling the TWOINTPTS function with the 
two cross-point intersections. The ISTANPT function can be used to determine 
which point is a tangent; the remaining two intersection points are then passed to 
TWOINTPTS. This logic is implemented in the THREEINTPTS function, with 
pseudo-code in Fig. 16. 

There is only one possible case for four intersection points, shown in Fig. 9. The 
two ellipse curves must cross at all four of the intersection points, resulting in a 
partial overlap. The overlap area consists of two segments from each ellipse, and a 
central convex quadrilateral. For the approach presented here, the four intersection 
points are sorted ascending in a counter-clockwise order around the first ellipse. 




a:^midRT ^Xraid ' COs(932 - 'fl) + Vmid ' sin((/3i - ip2) + /i2TR 

ymidRT =a:^mid ' sin((^2 - y^i) + 2/mid • cos((pi - IP2) + fexR 
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Figure 9. When three intersection points exist, one must be a 
tangent, and the ehipse curves must cross at the other two points, 
always resuhing in a partial overlap. When four intersection points 
exist, the ellipse curves must cross at all four points, again resulting 
in a partial overlap 



The ordered set of intersection points is (xi, j/i), (x2, 1/2), (2:3, ys) and {x4, yn). 
The ordering allows a direct calculation of the quadrilateral area. The standard 
formula uses the cross-product of the two diagonals: 

area =i |(a;3 - xi.yz - yi) x [x^ - 2:2,2/4 - ^2)! 

2 (40) 
= 2 1(^3 - 2:1) • (2/4 - 2/2) - (2:4 - 2:2) • (2:3 - 2:1)1 

The point ordering also simplifies the search for the appropriate segments of each 
ellipse that contribute to the overlap area. 

Suppose that the first two sorted points (2:1, 2/1) and (2:2, 2/2) demarcate a segment 
of the first ellipse that contributes to the overlap area, as shown in Fig. 9 and Fig. 10. 
ft follows that the contributing segments from the first ellipse are between (xi, 2/1) 
and (2:2, 2/2)7 and also between (0:3, 2/3) and (2:4, 2/4). In this case, the contributing 
segments from the second ellipse are between (0:2, 2/2) and (2:3, 2/3), and between 
(2:4, 2/4) and (2:1, 2/1). To determine which segments contribute to the overlap area, 
it suffices to test whether a point midway between (xi, 2/1) and (2:2, 2/2) is inside 
or outside the second ellipse. The segment algorithm is used for each of the four 
areas, and added to the quadrilateral to obtain the total overlap area. 




Figure 10. Overlap Area with four intersection points (Case 4-1). 
The overlap area consists of two segments from each ellipse, and a 
central convex quadrilateral. 
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An implementation of the ELLIPSE_ELLIPSE_OVERLAP algorithm in c-code 
is shown in Appendix6. The code compiles under Cygwin-1.7.7-1, and returns the 
following values for the test cases presented above in Fig. 6, Fig. 7 , Fig. 8 and 
Fig. 9: 



Listing 5. Return values for the test cases presented above in 
Fig. 6, Fig. 7 , Fig. 8 and Fig. 9. 



138 cc call\_cc.c c 11 i p s c \ _ c 1 1 i p s c \ _o v c r 1 a p . c — o i- a 1 I \ _cc . cxe 
139 

140 . / c al 1 \ _co 

141 

142 Calling c 1 1 i p s c \ _ c 1 1 i p s c \ _o v e r lap . c 

143 
144 
145 

146 Case 0-1: area = 6.28318531, rcturn.valuo = 111 

147 

148 ellipse 2 area by pi*a2*b2 = 6.28318531 

149 

150 Case 0-2: area = 6.28318531, r e t u r n_ v a 1 uo = 110 

151 

152 ellipse 1 area by pi*al*bl = 6.28318531 

153 

154 Case 0-3: area = 0.00000000, return.valuc = 103 

155 

156 Ellipses are disjoint , ovelap area — 0.0 

157 

158 

159 

160 Case 1-1: area = 6.28318531, return.valuc = 111 

161 

162 ellipse 2 area by pi*a2*b2 = 6.28318531 

163 

164 Case 1-2: area = 6.28318531, return.valuc = 110 

165 

166 ellipse 1 area by pi*albl = 6.28318531 

167 

168 Case 1-3: area = -0.00000000, return.valuc = 107 

169 

170 Ellipses are disjoint, ovelap area — 0.0 

171 

172 

173 

174 Case 2-1: area = 10.60055478, return.valuc = 109 

175 

176 ellipse 2 area by pi*a2*b2 = 10.60287521 

177 

178 Case 2-2: area = 6.28318531, return.valuc = 110 

179 

180 ellipse 1 area by pi*albl = 6.28318531 

181 

182 Case 2-3: area = 3.82254574, return.valuc = 107 

183 

184 

185 

186 Case 3-1: area = 7.55370392, return.valuc = 107 

187 

188 Case 3-2: area = 5.67996234, return.valuc = 107 

189 

190 

191 

192 Case 4-1: area = 16.93791852, return.valuc = 109 
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Listing 6. The ELLIPSEJ]LLIPSE_OVERLAP algorithm is 
shown for calculating the overlap area between two general el- 
lipses. The algorithm calls several supporting functions, including 
the polynomial solvers BIQUADROOTS, CUBICROOTS and 
QUADROOTS, from CACM Algorithm 326 [2] . The remaining 
functions are outlined in figures below. 



193 

194 (Area, Code) <- ELLIPSE.ELLIPSE.OVERLAP ( Al , Bl , HI , Kl , 1 , A2 , B2 , H2 , K2 , i/j2) 
195 

196 do if (Al = or Bl = 0) OR (A2 = or B2 = 0) 
197 

198 then return (-1, ERROR J^LLIPSEJARAMETERS) : DATA CHECK 

199 

200 do if {\ipl\ > 2it) 
201 

202 then ipl -i— (cplmodulo 2it- ) 

203 

204 do if {\tp2\ > 2tv) 
205 

206 then i,32 <— ((^2modulo 27r ) 

207 

208 H2.TR ^ (H2 - Hl)*cos (tpl) + (K2 - Kl)*sin (ipl) iTRANSfROT ELL2 
209 

210 K2.TR 4- (HI - H2)*sin (i^I) + (K2 - Kl)*cos (^pl) 
211 

212 ip2B, <- ip2 ipl 

213 

214 do if (|ip2_R| > 27r) 
215 

216 tlien ip2R -i^ (<p2Rmodulo2jr) 

217 

218 AA <- cos^(i^2R)/A2^ + sin^ (i^2R) /B2^ :BUILD\, IMPLICIT\ , COEFFS ELL2TR 
219 

220 BB <~ 2*cos (i^2R)*sin {tfi2R)/A2^ 2*cos ((,32R)*sin (ip2R)/B2^ 

221 

222 CC i- sin^ (i^2R) /A2^ + cos^ (ip2R) /B2^ 
223 

224 DD ^ -2*cos (ip2R)*(cos (¥>2R)*H2.TR + sin ( ip2R) *K2.TR) /A2^ 
225 

226 - 2*sin (¥J2R)*(sin ( i^2R) *H2.TR - cos ( i/j2R) *K2.TR ) /B2^ 

227 

228 EE ^ -2* sin (ip2R)* (cos (¥>2R)*H2.TR + sin ( ip2R) *K2.TR) /A2^ 
229 

230 + 2*cos (¥>2R)* (sin ( ip2R) * H2.TR - cos ( ip2R) *K2.TR) /B2^ 

231 

232 FF 4- (-cos (ip2R)*H2.TR - sin (ip2R) *K2.TR) ^ /A2^ 
233 

234 + (sin ((p2R) *H2.TR - cos ( ¥'2R) *K2.TR) ^ /B2^ - 1 

235 

236 : BUILD QUARTIC POLYNOMIAL COEFFICIENTS FROM THE TWO ELLIPSE EQNS 

237 

238 cy[4] ^ Al^+AA^ + Bl^ * ( Al^ * (BB^ - 2*AA*CC)+ B1^*CC^) 
239 

240 cy[3] ^ 2*BI*(B1^*CC*EE + A1^*(BB*DD - AA*EE) ) 
241 

242 cy[2] ^ AI^ * ( ( Bl^ * ( 2 * AA*CC — BB^ ) + DD^ - 2*AA*FF) 
243 

244 -2*AI^*AA^ + B1^*(2*CC*FF + EE^ ) 

245 

246 cy[l] ^ 2*BI*(AI^*(AA*EE — BB*DD) + EE*FF) 
247 

248 cy[0] i- (A1*(A1*AA — ^DD) + FF) * ( Al * ( A1*AA + DD) + FF) 
249 

250 py[0] ^ 1 
251 

252 do if (|cy[4]|> 0) : SOLVE QUARTIC EQ 

253 

254 then for i <- l.i 3 by 1 
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255 






256 




py[4- i ] <- cy [ i ] / cy [4] 


257 






258 




r [ ] [ ] <- BIQUADROOTS ( py [ ] ) 


259 






260 




n roots 4— 4 


261 






262 


else 


if (|cy[3]|>0) 


263 






264 


then 


for i 4- I o 2 by 1 


265 






266 




Py[3- i ] ^ cy [ i ] / cy [3] 


267 






268 




r [ ] [ ] <- CUBICROOTS ( py [ ] ) 


269 






270 




n roots -f- 3 


271 






272 


else 


if (|ci;[2]|>0) 


273 






274 


then 


for i <— 1 by 1 


275 






276 




py[2- i ] ^ cy [ i ]/cy [2] 


277 






278 




r [ ] [ ] ^ QUADROOTS ( py [ ] ) 


279 






280 




nroots 2 


281 






282 


else 


ll (|ci;[l]|>0) 


283 






284 


then 


r[l][l] +- (-cy[0]/cy [1]) 


285 






286 




r[2][l] ^ 


287 






288 




nroots 1 


289 






290 


else 




291 






292 




nroots 


293 






294 


nychk 





295 






296 


for i 


1 (o nroots by 1 


297 






298 


do if 


(|r[2]H| < EPS) 


299 






300 


then nychk ^ nychk + 1 


301 






302 




ychk[ nychk] i— r[l][i]* 


303 






304 


for j <- 


2 ( o nychk by 1 


305 






306 


tnipO 


ychk [ j ] 


307 






308 


for k 


^ (j — 1) to 1 by -1 


309 






310 


do 


if (ychk [k] ^ tmpO) 


311 






312 




then break 


313 






314 




else ychk[k + l] ychk[k] 


315 






316 


ychk[k + l] <- tmpO 


317 






318 


n i n t p t s 


^ 


319 






320 


for i 


1 1 o nychk by 1 


321 






322 


do if 


( ( 2 > 1 ) and ( |yc/iA;[i] — ychk['i 


323 






324 


then continue 


325 






326 


do if 


i\vchkli\\ > -Bl) 


327 







: SOLVE CUBIC EQ 



: SOLVE QUADRATIC EQ 



: SOLVE LINEAR EQ 



: COMPLETELY DEGENERATE EQ 



: IDENTIFY REAL ROOTS 



:SORT REAL ROOTS 



:FIND INTERSECTION POINTS 
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328 then xl •<— 

329 

330 else xl<-? Al*sqrt (1.0 - ychk [ i ] ^ /Bl^ ) 

331 

332 x2 <- -xl 

333 

334 do if {\elUpse2trixl,ychkli],AA,BB,CC,DD,EE,FF)\< EPS/2) 

335 

336 then nintpts -f- nintpts + 1 

337 

338 do if (nintpts > 4) 

339 

340 then return (-1, ERROR JNTERSECTIONJTS ) 

341 

342 xint [nintpts] -f— xl 

343 

344 yint[nintpts] -f— ychk[i] 

345 

346 do if ({\elUpse2tr{x2,ychk[i], AA,BB,CC,DD, EE, FF)\ < EPS/2) 

347 

348 and (|£c2-a;l|> EPS/2)) 

349 

350 then nintpts -f— nintpts + 1 

351 

352 do if (nintpts > 4) 

353 

354 then return (-1, ERRORJNTERSECTION JTS ) 

355 

356 xint [nintpts] -ir- xl 

357 

358 y i n t [ n i n t p t s ] -i— yehk [ i ] 

359 

360 s wit eh (nintpts) : HANDLE ALL CASES FOR \# OF INTERSECTION PTS 

361 

362 case 0: 

363 

364 case 1: 

365 

366 (OvorlapArea ,Codo) 4- NOINTPTS ( Al , Bl , A2 , B2 , HI , Kl , H2 , K2 , AA, 

367 

368 BB,CC,DD,EE,FF) 
369 

370 return ( O verlapArea , Code ) 

371 

372 case 2: 

373 

374 Code <- istanpt ( x i n t [ 1 ] , y i n t [ 1 ] , Al , Bl , AA, BB , CC , DD, EE , FF) 

375 

376 do if (Code = TANGENTJOINT) 

377 

378 then ( O verlapArea , Code ) <- NOINTPTS ( Al , Bl , A2 , B2 , HI , Kl , 

379 

380 H2,K2,AA,BB,CC,DD,EE,FF) 
381 

382 else ( OverlapArea , Code) ^TWOINTPTS ( x i n t [ ] , y i n t [ ] , Al , 

383 

384 PHI.l ,A2,B2,H2.TR,K2.TR, PHI.2 , AA, BB, CC.DD, EE, FF) 

385 

386 return (OverlapArea , Code ) 

387 

388 case 3: 

389 

390 (OverlapArea , Code) ^THREEINTPTS ( xint , yint , Al , Bl , PHI.l , 

391 

392 A2,B2,H2.TR,K2.TR,PHI.2 , AA, BB , CC,DD, EE, FF) 

393 

394 return (OverlapArea , Code) 

395 

396 case 4: 

397 

398 (OverlapArea , Code) <- FOURINTPTS ( xint , yint , Al , Bl , PHI.l , 

399 

400 A2, B2,H2.TR,K2.TR,PHI.2 ,AA,BB,CC,DD,EE,FF) 



401 
402 
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return { O vcrlapArca , Code ) 



Listing 7. The NOINTPTS subroutine. If there are either or 1 
intersection points, this function determines whether one elhpse is 
contained within the other (Cases 0-1, 0-2, 1-1 and 1-2), or if the 
eUipses are disjoint (Cases 0-3 and 1-3). The function returns the 
appropriate overlap area, and a code describing which case was 
encountered. 



403 

404 (OvcrlapArca ,Codc) <- NOINTPTS ( Al , Bl , A2 , B2 , HI , Kl , H2.TR , K2.TR , AA, 
405 

406 BB,CC,DD,EE,FF) 
407 

408 rolsizc <- A1*B1 - A2*B2 
409 

410 do if (rclsizc > 0) 
411 

412 then do if (((H2.TR*H2.TR)/(A1*A1)+(K2.TR*K2.TR)/(B1*B1) ) < 1.0) 

413 

414 then return (7r*A2*B2 , ELLIPSE2.INSIDE.ELLIPSE1 ) 

415 

416 else return (0, DIS JOINT.ELLIPSES ) 

417 

418 else do if (relsizc < 0) 

419 

420 then do if (FF < 0) 

421 

422 then return ( 7r * Al* Bl , ELLIPSE1.INSIDE.ELLIPSE2 ) 

423 

424 else return (0, DISJOINT.ELLIPSES ) 

425 

426 else do if ((HI = H2.TR) AND (Kl = K2.TR) ) 

427 

428 then return (7r*Al*Bl, ELLIPSES J\.REJDENTICAL ) 

429 

430 else return (-1, ERROR-CALCULATIONS 

Listing 8. The TWOINTPTS subroutine. If there are 2 intersec- 
tion points where the ehipse curves cross (Case 2-3), this function 
uses the elhpse sector algorithm to determine the contribution of 
each ellipse to the total overlap area. The function returns the 
appropriate overlap area, and a code indicating two intersection 
points. 

431 (OverlapArea ,Code) <- TWOINTPTS ( x i n t [ ] , y i n t [ ] , Al , Bl , ¥J 1 , A2 , B2 , H2.TR , 
432 

433 K2.TR,ip2 ,AA,BB,CC,DD,EE,FF) 

434 

435 do if (|a;[l]|> Al) lAVOID INVERSE TRIG ERRORS 

436 

437 then do if (x[l] < 0) 

438 

439 then x [ 1 ] <- -Al 

440 

441 els c X [ 1 ] •(- Al 

442 

443 do if (y[l] < 0) : FIND PARAMETRIC ANGLE FOR (x[l], y[l]) 

444 

445 tlicn ei <- 2iT arccos (x[l]/Al) 

446 

447 else ei arccos (x[l]/Al) 

448 

449 do if (|x[2]|> Al) : AVOID INVERSE TRIG ERRORS 
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450 

451 then do if (x[2] < 0) 

452 

453 then x [ 2 ] •(- -Al 

454 

455 e 1 s c X [ 2 ] <- Al 

456 

457 do if (y[2] < 0) : FIND PARAMETRIC ANGLE FOR (x[2], y[2]) 

458 

459 Uion 92 ^ Itt arccos (x[2]/Al) 

460 

461 else 62 ^ arecos (x[2]/Al) 

462 

463 do if (e\ > 62) :GO CCW FROM 61 TO\ , 92 

464 

465 tlien tmp <- 9 1 , S 1 <- S 2 , 92 <- tmp 

466 

467 xmid <- Al*cos {(61 + 62) /2} 
468 

469 ymid i- Bl*sin ((91 + 62) /2} 
470 

471 do if (AA*xinid^+BB*xmid*ymid+CC*yniid^4DD*xmid+EE*ymid+FF > 0) 
472 

473 then tmp <- 91, 91-1- 62, 92 tmp 

474 

475 do if (91 > 62) : SEGMENT ALGORITHM FOR ELLIPSE 1 

476 

477 Llien 61 1 61 ~ 27r 

478 

479 do if ((62 ~ 91) > rr) 
480 

481 til en trsign 1 

482 

483 else trsign — 1 

484 

485 areal ■(- (Al*Bl*(e2 - 91) + t r s i g n * \ t ox t b ar x[l]*y[2] - x [ 2 ] * y [ 1 ] ) \ toxtbar 
/2 

486 

487 xl\.tr <- (x[l] - H2\.TR) * cos (i^l — ip2) + (y[l] - K2\.TR) * sin ((p2 — ipl) 
488 

489 yl\.tr <- (x[l] - H2\.TR) * sin (c^l — ip2) + (y[l] - K2\.TR) * cos ((pi — ifi2) 
490 

491 x2\.tr ^ (x[2] - H2\.TR) * c os ( ip 1 ip2) + (y [2] - K2\.TR) * sin ((p2 ifil) 

492 

493 y2\.tr ? (x[2] - H2\.TR) * s i n ( y 1 — ip2) + (y [2] - K2\.TR) * cos ((^1 — ip2) 
494 

495 do if (|a;l.tr|> A2) : AVOID INVERSE TRIG ERRORS 

496 

497 then do if (xl.tr < 0) 

498 

499 then xl.tr -A2 

500 

501 else xl.tr ^ A2 

502 

503 do if (yl.tr < 0) : FIND PARAMETRIC ANGLE FOR (xl.tr, yl.tr) 

504 

505 then 91 <- 27r arccos (xl.tr/A2) 

506 

507 else 91 <- arccos (xl.tr/A2) 

508 

509 do if (\x2tr\ > A2) : AVOID INVERSE TRIG ERRORS 

510 

511 then do if (x2.tr < 0) 

512 

513 then x2.tr i- -A2 

514 

515 else x2.tr i- A2 

516 

517 do if (y2.tr < 0) : FIND PARAMETRIC ANGLE FOR (x2.tr, y2.tr) 

518 

519 then 92 ^ 27r arccos (x2.tr/A2) 

520 

521 else 92 ^ arccos (x2.tr/A2) 
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522 

523 do if {01 > 02) :GO CCW FROM 01 TO\ , 02 

524 

525 then tmp •(- S 1 , S 1 <- S 2 , 612 <- tmp 

526 

527 xmid <- A2*cos {{01 + e2)/2) 
528 

529 ymid <- B2*sin {{01 + 02)12) 
530 

531 xmid_rt — xmid* c os ( :p2 1 ) + ymid * s i n ( c/; 1 ip2) + H2_TR 

532 

533 ymid_rt — xmid*sin(c^2 1 ) + ymid*cos(:^2 ^1) + K2_TR 

534 

535 do if ( xmid_rt^/Al^ + ymid_rt^/Bl^ > 1) 
536 

537 Uicu tmp <- ei, ei •(- 62, 612 <- tmp 

538 

539 do if {01 > 02) : SEGMENT ALGORTTHM FOR ELLIPSE 2 

540 

541 tlicii ei •(- 91 - 2ir 

542 

543 do if {{02 - 01) > ir) 
544 

545 t ii c n t r s i g n 1 

546 

547 else trsign — 1 

548 

549 area2 <- (A2*B2*(6/2 - 01) 
550 

551 + tr si gn ^\xltr ^ y2tr — x2^tr ^ yltr)\ /2 
552 

553 ret 111 11 (areal + area2 , TWO JNTERSECTION.POINTS ) 



Listing 9. The THREEINTPTS subroutine. When there are 
three intersection points, one of the points must be a tangent 
point, and the eUipses must cross at the other two points. For 
the purpose of determining overlap area, the TWOINTPTS 
function can be used with the two cross-point intersections. The 
ISTANPT function can be used to determine which point is a 
tangent; the remaining two intersection points are then passed to 
TWOINTPTS. The function returns the appropriate overlap area, 
and a code indicating three intersection points. 



554 OverlapArea , Code) ^THREEINTPTS ( x i n t [ ] , y i n t [ ] , Al , Bl , 1 , A2 , B2 , H2.TR , 
555 

556 K2.TR,ip2 ,AA,BB,CC,DD,EE,FF) 

557 tanpts 
558 

559 for i 1 lo nych]< by 1 

560 

561 eode ^ ISTANPT ISTANPT ( x [ i ] , y [ i ] , Al , Bl , AA, BB, CC,DD, EE, FF) 
562 

563 do if (code = TANGENT JOINT) 
564 

565 then tanpts tanpts + 1 

566 

567 tanindex -ir- i 

568 

569 do if NOT (tanpts = 1) 
570 

571 then return (-1, ERRORJNTERSECTION J'OINTS ) 
572 

573 switch (tanindex) : STORE THE INTERSECTION POINTS 

574 

575 case 1: : TANGENT POINT IS IN (x[l], y[l]) 

576 

577 xint[l] xint[3] 
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578 

579 yint[l]-(-yint[3] 
580 

581 case 2: : TANGENT POINT IS IN (x[2], y[2]) 

582 

583 xint[2]-(-xint[3] 
584 

585 y i n t [ 2 ] y i n t [ 3 ] 

586 

587 (OvorlapAroa , code) <-TWOINTPTS ( x i n t [ ] , y i n t [ ] , Al , Bl , ¥J 1 , A2 , B2 , H2.TR , 
588 

589 K2.TR,ip2 ,AA,BB,CC,DD,EE,FF) 

590 

591 letuin ( Overlap Area , THREEJNTERSECTION JOINTS) 



LISTING 10. The FOURINTPTS subroutine. When there are four 
intersection points, the eUipse curves must cross at ah four points. 
A partial overlap area exists, consisting of two segments from 
each ellipse and a central quadrilateral. The function returns the 
appropriate overlap area, and a code indicating four intersection 
points. 



592 


verlapArea , Code) <- FOURINTPTS (xint[] , 


yint [] ,A1,B1,¥j1 ,A2,B2,H2.TR, 


593 






594 


K2.TR,^2 ,AA,BB,CC,DD,EE,FF) 


595 






596 


for i <— 1 I o 4 by 1 


: AVOID INVERSE TRIG ERRORS 


597 






598 


do if (|xint[i]| > Al) 




599 






600 


then do if ( xint [ i ] < 0) 




601 






602 


tlien xint[i] -i— — Al 




603 






604 


else xint[i] -t— Al 




605 






606 


do if (yint [ i ] < 0) 


:FIND PARAMETRIC ANGLES 


607 






608 


then ^[i] ^ 2-jr arccos (xint[ 


i]/Al) 


609 






610 


else ^[i] arccos (xint[i]/Al) 




611 






612 


for j •<- 2 to 4 by 1 


:PUT POINTS IN CCW ORDER 


613 






614 


tmpO ^ e[j ] 




615 






616 


tinpl xint [ j ] 




617 






618 


tmp2 <- yint [ j ] 




619 






620 


for k <- ( j -1) to 1 by -1 


: INSERTION SORT BY ANGLE 


621 






622 


do if (e[k] <= tmpO) 




623 






624 


then break 




625 






626 


else e[k+l] ^ elk] 




627 






628 


xint [k + 1] <- xint [k] 




629 






630 


yint [k + 1] 4- yint [k] 




631 






632 


a real ^ (|(xint[3] xint [1]) *(yint [ 


4] — yint [2]) — 


633 






634 


xint [4] - xint [2] )*( yint [3] yint[l]) 


1 /2) :QUAD AREA 


635 






636 


for i 1 f o 4 by 1 


:TRANSLATEfROTATE ELLIPSE 2 


637 






638 


xint_tr[i] <- (xint[i] H2_TR)*cos (i^l ip2) 
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639 

640 + (yint[i] — K2.TR)*sin {ip2 — (pi) 

641 

642 yint_tr[i] <- (xint[i] H2_TR)*sin (ipl ip2) 

643 

644 + (yint[i] K2.TR)*cos {ipl ip2) 

645 

646 do if {\xinttr[i]\ > A2) : AVOID INVERSE TRIG ERRORS 

647 

648 then do if (xint_tr[i] < 0) 

649 

650 thcuxint_tr[i]<— — A2 

651 

652 else xint_tr[i] <- A2 

653 

654 do if (yint.tr [i] < 0) : FIND PARAM ANGLES FOR (xint.tr, yint.tr) 

655 

656 then e.tr[i] <— 2tv arccos ( x i n t .t r [ i ] / A2 ) 

657 

658 else e.tr[i] <- arceos ( x i n t . t r [ i ] / A2 ) 

659 

660 xmid •(- Al*cos ((61 + 02) /2) 
661 

662 yniid 4- Bl*sin ((91 + e2)/2) 
663 

664 do if (AA*xmid^+BB*xmid*ymid+CC*yniid^4DD*xmid+EE*ymid+FF < 0) 
665 

666 then area2 = (A1*B1*(0[2] - [1\) 

667 

668 - I ( xint [1] * yint [2] - x i n t [ 2 ] * y i n t [ 1 ] ) | ) /2 
669 

670 areas = (A1*B1*(0[4] - 0[3]) 

671 

672 - I ( xint [3] * yint [4] - x i n t [ 4 ] * y i n t [ 3 ] ) | ) /2 
673 

674 area4 = ( A2* B2 * ( 9 . t r [ 3 ] - e.tr[2]) 

675 

676 - I ( xint\ .tr [2] * yint.tr [3] - x i n t . t r [ 3 ] * y i n t . t r [ 2 ] ) | )/2 
677 

678 areas = ( A2* B2 * ( S . t r [ 1 ] - e.tr[4] - twopi)) 

679 

680 - I ( xint.tr [4] * yint.tr [1] - x i n t .t r [ 1 ] * y i n t .t r [ 4 ] ) | / 2) 
681 

682 else area2 = (A1*B1*(9[3] - 9 [2]) 

683 

684 - I ( xint [2] * yint [3] - x i n t [ 3 ] * y i n t [ 2 ] ) | ) /2 

685 

686 areas = (A1*B1*(0[1] - (0[4] - twopi)) 

687 

688 - I ( xint [4] * yint [ 1 ] - x i n t [ 1 ] * y i n t [ 4 ] ) | ) /2 
689 

690 area4 = ( A2* B2 * ( . t r [ 2 ] - e.tr[l]) 

691 

692 - I ( xint.tr [1] * yint.tr [2] - x i n t . t r [ 2 ] * y i n t . t r [ 1 ] ) | ) /2 
693 

694 area5 = ( A2* B2 * ( 6* . t r [ 4 ] - e.tr[3]) 

695 

696 - I ( xint.tr [3] * yint.tr [4] - x i n t . t r [ 4 ] * y i n t . t r [ 3 ] ) | ) /2 
697 

698 return ( areal + aroa2 + area3+area4+aroa5 , FOURJNTERSECTION JOINTS ) 

Listing 11. The ISTANPT subroutine. Given an intersection 
point (a;, y) that satisfies both EUipse Eq.21 and Elhpse Eq. 22, 
the function determines whether the two elhpse curves are tangent 
at (i, y), or if the ellipse curves cross at (x, y). 

699 Code ^ ISTANPT (x , y , Al , Bl , AA, BB, CC,DD,EE, FF) 
700 

701 do if (|a;| > Al) : AVOID INVERSE TRIG ERRORS 

702 
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703 then do if x < 

704 

705 then X ^ -Al 

706 

707 else x i- Al 

708 

709 do if (y < 0) : FIND PARAMETRIC ANGLE FOR (x, y) 

710 

711 then 6 <— 27r arccos (x/Al) 

712 

713 else 6 <— arccos (x/Al) 

714 

715 branch ^ v(x^ + ) : DETERMINE PERTURBATION ANGLE 

716 

717 .lo if (branch < 100*EPS) 
718 

719 then eps.radian <- 2*EPS 

720 

721 else eps_radian i~ arcsin ( 2 * EPS/ branch ) 

722 

723 xl <- Al*cos (e + eps.radian) iCREATE TEST POINTS ON EACH SIDE 
724 

725 yl ^ Bl*cos (e + eps.radian) :OF THE INPUT POINT (x, y) 

726 

727 x2 Al*cos {0 — eps.radian) 
728 

729 y2 -f— Bl*cos {0 — eps.radian) 
730 

731 tostl <- AA*xl^+BB*xl*yl+CC*yl^-ff)D*xl+EE*yl+FF 
732 

733 tost2 <- AA*x2^+BB*x2*y2+CC*y2^-H)D*x2+EE*y2+FF 
734 

735 do if ( test 1* test 2 > 0) 
736 

737 then return TANGENT JOINT 

738 

739 else return INTERSECTION JOINT 



Listing 12. C-SOURCE CODE FOR ELLIPSE.SEGMENT 

4. APPENDIX A. 

740 
741 

742 /* 
743 

744 * 
745 

746 * Function: double elli'p s e e gment 
747 

748 * 
749 

750 * Purpose : Given the parameters of an ellipse and two points that lie on 
751 

752 * the ellipse , this function calculates the ellipse segment 

area 

753 

754 * between the secant line and the ellipse . Points are input as 

755 

756 * (XI, yi) and (X2, Y2) , and the segment area is defined to be 

757 

758 * between the secant line and the ellipse from the first point 

759 

760 * (XI, Yl) to the second point (X2, Y2) in the counter- 

clockwise 

761 

762 * direction. 
763 

764 * 
765 
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766 
767 
768 
769 
770 
771 
772 
773 
774 
775 
776 
777 
778 
779 
780 
781 
782 
783 
784 
785 
786 
787 
788 
789 
790 
791 
792 
793 
794 
795 
796 
797 
798 
799 
800 
801 
802 

803 
804 
805 
806 



^ Reference : Hughes and Chraibi (2011) . Calculating Ellipse Overlap Areas 

* Dependencies : math . h for calls to trig and absolute value functions 

* pro gram^constants . h error message codes and constants 



Inputs . 



Outputs : 



Return : 



1 . 


double 


A 


e 


llip s e 


serni— axis length in 


X— direction 


2. 


double 


B 


e 


Hip s e 


seini— axis length in 


y— direction 


3. 


double 


XI 


X' 


-value 


of the first point 


on the ellipse 


4- 


double 


Yl 


y- 


-value 


of the first point 


on the ellipse 


5. 


double 


X2 


X- 


-value 


of the second p oint 


on the ellipse 


6. 


double 


Y2 


y- 


-value 


of the second p oint 


on the ellipse 



1 . int -^MessageCode stores diagnostic inform ation 

integ er codes in program^constants . h 

The V alue of the ellipse segment area : 

— 1.0 is returned in case of an error with input data 



// 



808 //== INCLUDE ANSI C SYSTEM AND USER^DEFINED HEADER FILES 



809 

810 // 



811 
812 
813 
814 
815 
816 

817 
818 
819 
820 
821 
822 
823 
824 
825 
826 
827 
828 

829 
830 



^^include " p r ogr am_c o ns t an t s . h" 

double ellipsc_segmcnt (double A, double B, double XI, double Yl , double X2 

double Y2 , inl * McssagcCode ) 



{ 



double thetal; / / — p arametric angle of the first p oint 
double theta2; / / — p arametric angle of the second point 
double trsign; / / — sign of the triangle area 

double pi — 2.0 * asin \eqref{GrindEQ 1_0_}; / / a maximumr- 

pr e c i s i o n value of p i 



double twopi — 2.0 * pi ; 



/ / — a maximuiTb- p r ecision v alue of 2^ pi 
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831 
832 
833 

834 // Check the data first 

835 

836 // Each of the ellipse axis lengths must he positive 

837 

838 if (!(A > 0.0) \tcxtbar \tcxtbar !(B > 0.0)) 

839 

840 { 
841 

842 (* MossagoCodc ) = ERROR^LLIPSE PARAMETERS ; 

843 

844 return -1.0; 

845 

846 } 

847 

848 

849 

850 // Points must be on the ellipse , within EPS, which is defined 

851 

852 // in the header file program^constants.h 

853 

854 if ( (fabs ((X1*X1) /(A*A) + ( Y1*Y1 ) / (B*B) - 1.0) > EPS) toxtbar 

t e X t b a r 

855 

856 (fabs ((X2*X2) /(A*A) + (Y2*Y2) / (B*B) - 1.0) > EPS) ) 

857 

858 { 
859 

860 (*McssagcCodG) = ERROR JOINTS J^OT.ON^LLIPSE ; 

861 

862 return -1.0; 

863 

864 } 

865 

866 

867 

868 // Avoid inverse trig calculation errors: there could be an error 

869 

870 // if \textbar Xl/A\textbar > 1.0 or \textbar X2/A\texthar > 1.0 

when calling acos() 

871 

872 // // execution arrives here, then the point is on the ellipse 

873 

874 // within EPS. Try to adjust the value of XI or X2 before giving 

875 

876 // — up on the area calculation 
877 

878 if (fabs (Xl)/A > 1.0) 

879 

880 { 
881 

882 // — if execution arrives here, already know that \textbar Xl\ 
textbar > A 

883 

884 if ((fabs (XI) - A) > EPS) 

885 

886 { 
887 

888 / / — i f XI is not close to A or —A . then give up 



890 (* MessageCodo ) = ERROR\.INVERSE\.TRIG ; 

891 

892 return -1.0; 

893 

894 } 
895 

896 else 
897 

898 { 
899 

900 // — nudge XI back to A or —A, so acos() will work 
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901 

902 XI = (XI < 0) ? -A : A; 

903 

904 } 
905 

906 } 

907 

908 

909 

910 if (fabs (X2) /A > 1.0) 

911 

912 { 

913 

914 // — if execution arrives here, already know that \texthar X2\ 
textb ar > A 

915 

916 if ((fabs (X2) - A) > EPS) 

917 

918 { 
919 

920 / / — i f X2 is not close to A or —A . then give up 

921 

922 (»MessageCode) = ERRORJNVERSE.TRIG ; 

923 

924 roturn— 1.0; 

925 

926 } 
927 

928 else 
929 

930 { 
931 

932 // — nudge X2 back to A or ~A, so acos() will work 
933 

934 X2 = (X2 < 0) ? -A : A; 

935 

936 } 
937 

938 } 

939 

940 

941 

942 // Calculate the parametric angles on the ellipse 

943 

944 // The parametric angles depend on the quadrant where each point 

945 

946 // is located. See Table 1 in the reference. 

947 

948 if (Yl < 0.0) // — Quadrant III or IV 

949 

950 thctal = twopi — acos (XI / A) ; 

951 

952 else // — Quadrant I or II 

953 

954 thctal = acos (XI / A) ; 

955 

956 

957 

958 if (Y2 < 0.0) // — Quadrant III or IV 

959 

960 thcta2 = twopi — acos (X2 / A) ; 

961 

962 else // — Quadrant I or II 

963 

964 thcta2 = acos (X2 / A) ; 

965 

966 

967 

968 // need to start the algorithm with thctal < theta2 

969 

970 if (thctal > theta2) 

971 

972 thctal — — twopi; 
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973 




974 




975 




976 


// if the integration angle is less than pi. subtract the triangle 


977 




978 


// area from the sector , otherwise add the triangle area. 


979 




980 


if ((thcta2 — thetal) > pi) 


981 




982 


trsign — 1.0; 


983 




984 


else 


985 




986 


trsign — — 1.0; 


987 




988 
989 




990 


// The ellipse segment is the area between the line and the ellipse , 


991 




992 


// calculated by finding the area of the radial sector minus the 




area 


993 




994 


// of the triangle created by the center of the ellipse and the two 


995 




996 


// points. First term is for the ellipse sector; second term is for 


997 




998 


// the triangle between the points and the origin. Area calculation 


999 




1000 


// is described in the reference. 


1001 




1002 


(*MossagoCodc) = NORMAL.TERMINATION : 


1003 




1004 


return ( . 5 * (A*B * ( t h e t a2 - thetal) + trsign*fabs (X1*Y2 - X2*Y1)) ); 


1005 




1006 




1007 




1008 } 





Listing 13. C-SOURCE 
LIPSE.LINE.OVERLAP 



CODE 



FOR 



EL- 



1010 
1011 
1012 



1013 
1014 
1015 
1016 
1017 
1018 
1019 
1020 
1021 
1022 
1023 
1024 

1025 
1026 
1027 
1028 
1029 
1030 
1031 
1032 
1033 
1034 



5. APPENDIX B. 



Function: double ellip s e ^lin e ^ o v erl ap 



Purpose : 



tangent 



Y2) 



Given the parameters of an ellipse and two points on a line , 
this function calculates the area between the two curves. If 
the line does not cross the ellipse, or if the line is 

to the ellipse , then this function returns an area of 0.0 
If the line intersects the ellipse at two points , then the 
function returns the area between the secant line and the 
ellipse . The line is considered to have a direction from 
the first given p oint (XI . Yl ) to the second given p oint ( X2 , 
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1035 
1036 
1037 
1038 
1039 
1040 
1041 
1042 
1043 
1044 
1045 
1046 
1047 
1048 
1049 
1060 
1061 
1062 
1053 
1054 
1055 
1066 

1067 
1068 
1069 
1060 
1061 
1062 
1063 
1064 
1066 
1066 
1067 
1068 
1069 
1070 
1071 
1072 
1073 
1074 
1076 
1076 
1077 
1078 
1079 
1080 
1081 
1082 
1083 
1084 
1086 
1086 
1087 
1088 
1089 
1090 
1091 
1092 
1093 
1094 
1095 
1096 
1097 
1098 

1099 
1100 
1101 
1102 



This function determines where the line crosses the ellipse 
first , and where it crosses second. The area returned is 
between the secant line and the ellipse traversed counter- 
clockwise from the first intersection point to the second 
intersection point. 

Reference : Hughes and Chraibi (2011) , Calculating Ellipse Overlap Areas 

Dependencies : math . h for calls to trig and absolute value functions 
program-constants .h error message codes and constants 
ellipse-segment . c core algorithm for ellipse segment 



Inputs : 



1 . 


double 


PHI 


(JCW rotation angle of the ellipse , radians 


2. 


double 


A 


ellipse semi— axis length in x— direction 


3. 


double 


B 


ellipse semi— axis length in y— direction 


A- 


double 


H 


horizontal offset of ellipse center 


5. 


double 


K 


vertical offset of ellipse center 


6. 


double 


XI 


X— V alue of the first point on the line 


7. 


double 


Yl 


y— value of the first p oint on the line 


8. 


double 


X2 


X— value of the second point on the line 


9. 


double 


Y2 


y— value of the second point on the line 



Outputs : 



Return : 



1. int ^MessageCode returns diagnostic information 

integer codes in program-constants . h 

The v alue of the ellipse segment area : 

— 1.0 is returned in case of an error with the data or 
calculation 

0.0 is returned if the line does not cross the ellipse , or if 
the line is tangent to the ellipse 



*/ 



1103 
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1104 //== DEFINE PROGRAM CONSTANTS 



1113 
1114 



1115 
1116 



1105 

1106 // 



1108 9^include " p r ogr am_c ons t an t s . h" // error message codes and constants 

1109 
1110 
1111 

1112 // 



//== DEPENDENT FUNCTIONS 



1118 doLiljlc t ex t b f { c II i ps c _s e g m e n t } (double A, double B, double XI, double Yl , 
double X2, 

1119 

1120 double Y2 , int * MessageCode ) ; 

1121 

1122 

1123 

1124 double \ t e X t b f { e 1 1 i p s e _ 1 i n e _ o V c r 1 a p } (double PHI, double A, double B, 
double H , 

1125 
1126 
1127 
1128 
1129 

1130 \{ 
1131 

1132 // 



double K, double XI, double Yl , double X2 , 
double Y2 , int * MessageCode ) 



1133 
1134 



1135 
1136 



//== DEFINE LOCAL VARIABLES 
// 



1137 
1138 
1139 
1140 
1141 
1142 
1143 
1144 
1145 
1146 

1147 
1148 



double 


XIO; 


// — Translated , 


Rotated 


X— V alue of the 


first 1 


'J oint 


double 


YIO; 


// — Translated , 


Rotated 


y— V alue of the 


first 1 


'J oint 


d o u l3 1 e 


X20; 


// — - Translated , 


Rotated 


x~ V alue of the 


second 


p oint 


double 


Y20; 


// — ' Translated , 


Rotated 


y— V alue of the 


second 


p oint 


double cosphi 

mult i p I c 


^ textbf{cos} (PHI) 

c a I c s 


; //— 


store cos ( PHI ) 


to avoid 


double 


s i n p h i 


^ \ tcxtbf{sin} (PHI 


) ; /A- 


~ store sin ( PHI ) 


t o avo i d 



multiple c ales 



1150 


double 


m; 


//— 


line slope , 


calculated from input line slope 


1151 












1152 


double 


a , b , c ; 


//" 


~ quadratic 


equation co efficients a*3;\'{}2 -f b=¥x 


1153 


+ 


c 








1154 


double 


discrim ; 


//- 


— quadratic 


equationdis criminant — 4^ a^ c 


1155 












1156 


double 


xl , x2 ; 


/A 


-- X— values 


of intersection p oints 


1157 












1158 


double 


yl, y2; 


//- 


— y— values 


of intersection p oints 
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1159 
1160 
1161 
1162 
1163 
1164 
1165 
1166 
1167 
1168 
1169 
1170 
1171 
1172 
1173 
1174 
1175 
1176 
1177 
1178 
1179 
1180 
1181 
1182 
1183 
1184 
1185 
1186 
1187 
1188 
1189 
1190 
1191 
1192 

1193 
1194 
1195 
1196 
1197 
1198 
1199 
1200 
1201 
1202 
1203 
1204 
1205 
1206 
1207 
1208 
1209 
1210 
1211 
1212 
1213 
1214 
1215 
1216 
1217 
1218 

1219 
1220 
1221 
1222 
1223 
1224 
1225 
1226 
1227 
1228 
1229 



double mid_X ; // — midpoint of the rotated x— values on the line 

double thetalparm; // parametric angle of first point 

double theta2parm; // parametric angle of second point 

double xmidpoint; // x-— value midpoint of secant line 

double yinidpoint; // y— value midpoint of secant line 

double rootl , root2; // — temporary storage variables for roots 

double segmcnt.area ; // — stores the ellipse segment area 

// Check the data first 

// Each of the ellipse axis lengths must he positive 

if (!(A > 0.0) \tGxtbar \tGxtbar !(B > 0.0)) 
{ 

(♦Message Code) = ERROR JILLIPSEJARAMETERS ; 
return —1.0; 

} 



/ / The rotation angle for the ellipse should be b etween —2pi and 2pi 

(?) 

if ( (\textbf{fabs} (PHI) > (2.0*pi)) ) 
PHI ^ \ tcxtbf {fmod} {PHI, twopi); 



/ /— — For this numerical routine , the ellipse will he translated and 

/ / rotated so that it is centered at the origin and orient e d with 

// the coordinate axes . 

/ /— Then , the ellipse will have the implicit ( p olynoinial ) form of 
//— x\~{}2/A\-{}2 + y+2/B\-{}2 = 1 

// For the line . the given points are first translated by the amount 

/ / required to put the ellipse at the origin . e.g.. by (—H, —K) . 

// Then, the points are rotated by the amount required to orient 

/ / the ellipse with the coordinate axes . e.g., through the angle — 

PHI. 

XIO = cosphi*(Xl - H) + sinphi*(Yl - K) ; 
YIO = -sinphi*(Xl - H) + cosphi*(Yl - K) ; 
X20 = cosphi*(X2 - H) + sinphi*(Y2 - K) ; 
Y20 = -sinphi*(X2 - H) + cosphi*(Y2 - K); 
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1230 


//- 


— T^o deter7nin& if the line and ellipse intersect , solve the two 


1231 






^?oq 




— equdtiojis SI m uLtdneousLy , by substitutmg y — Y 1 (J -f- ( x — J(. 1 (J J 


1 






1234 




— and X — JilU -f- 7nxy^ ( y — Y 10 J tnt o the e 1 1 1 p s e equation , 


1235 






1236 




— which results m two quadratic equations m x . oee the reference 


1237 






1238 


//- 


— for derivations of the quadratic coefficients . 


1239 






1240 






1241 






1242 




If the 7iew line is not close to o etng vertical, then use the 


1243 






1244 


//- 


— first derivation 


1245 






1246 


1 I 










1 


r 

\ 




1249 






1250 












1252 




/ / ^ * ( 11 u^7n — 7n\ ^jii^-A.iuj * x 


1253 






1254 




( Y ltJ\ — ii^^nw Y H)^ .\1U -f- m\ ■[J-<;*Ai(/\ — -o\ xj'^J 


1255 






1256 




m ^ (Y20 - YIO) /{X20 - XIO) ; 


1257 






1258 




a — (B*B + A*A*m*m) /{A*A) ; 


1259 






1260 




2.0*(Y10*m - m*m*X10) ; 


1261 






1262 




c ^ (YIO^YIO - 2.0*m*Y10*X10 + m*m*X10*X10 - B*B) ; 


1263 






1264 


} 




1265 






1266 


//- 


— // the new line is do s e to b eing vertical, then use the 


1267 






1268 


//- 


— second derivation 


1269 






1270 


else if (\textbf{fabs} (Y20 - YIO) > EPS) 










{ 














//— ((A\'{}2 + B\'{}2*rn\-{}2)/(B\'{}2)) * y\-{}2 


1275 






1276 




// — 2* (X10*m - m\ '{} 2*Y10) * y 












// — (X10\~{}2 - 2*mY10*X10 + m\' {} 2*Y10\' {} 2 - A\-{}2) 


1279 






1280 




m = (X20 - XIO) /(Y20 - YIO) ; 


1281 






1282 




a = (A*A + B*B*m*m) /(B*B) ; 


1283 






1284 




b = 2.0*(X10*m- m*m*Y10) ; 


1285 






1286 




c = (X10*X10 - 2.0*m*Y10*X10 + m*in* YIO* YIO - A*A) ; 


1287 






1288 


} 




1289 






1290 


/A 


— // the two given points on the line are very close to g ether in 


1291 






1292 


//- 


— both X and y directions , then give up 


1293 






1294 


else 


1295 






1296 


{ 




1297 






1298 




(♦Message Code) = ERRORXINE JOINTS ; 


1299 






1300 




return —1.0; 


1301 






1302 


} 
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1303 






1304 






1305 






1306 


//- 


— Once the co efficients for the Quadratic Equation 


1307 






1308 


//- 


— known , the roots of the quadratic p olynornial will 


1309 






1310 


//- 


— the x— o r y— values of the p oints of intersection 


1311 






1312 


//- 


— and the ellipse. The dis criininant can be used to 


1313 






1314 


//- 


— which case has o ccurred for the given inputs : 


1315 






1316 


//- 


— 1 . dis cr <i 


1317 






1318 


//- 


— Quadratic has complex conjugate roots. 


1319 






1320 


//- 


— The line and ellipse do not intersect 


1321 






1322 


//- 


— 2. discr — 


1323 






1324 


//- 


— Quadratic has one repeated root 


1325 






1326 


//- 


— The line and ellipse intersect at only one 


1327 






1328 


/A 


— i.e.. the line is tang eiit to the ellipse 


1329 






1330 


//- 


— 3 . discr > 


1331 






1332 


//- 


— Quadratic has two distinct real roots 


1333 






1334 


//- 


— Th e line crosses the ellipse at two p oints 


1335 






1336 


discrim — b*b — 4.0*a*c; 


1337 






1338 


if 


(discrim < 0.0) 


1339 






1340 


{ 




1341 






1342 




/ / — Line and ellipse do not intersect 


1343 






1344 




(*Mcss age Code) ^ NO_INTERSECTION_POINTS ; 


1345 






1346 




return 0.0; 


1347 






1348 


} 




1349 






1350 


else if (discrim > 0.0) 


1351 






1352 


{ 




1353 






1354 




// — Two real roots exist , so calculate them 


1355 






1356 




/ / — The larger root is stored in root2 


1357 






1358 




rootl ^ (-b - \ tcxtbf { sqrt } (discrim)) / (2.0*a); 


1359 






1360 




root2 ^ (-b + \ tcxtbf { sqrt } (discrim)) / (2.0*a); 


1361 






1362 


} 




1363 






1364 


els 


c 


1365 






1366 


{ 




1367 






1368 




/ / — L i 71 e i s I a n g e n t t o the ell i p s e 


1369 






1370 




(*Mess age Code) ^ LINE_TANGENT_TO_ELLIPSE ; 


1371 






1372 




return 0.0; 


1373 






1374 


} 




1375 
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1376 
1377 

1378 // decide which roots go into which x or y values 

1379 

1380 if (\ textbf { fabs } (X20 - XIO) > EPS) // — roots are x-values 

1381 

1382 { 
1383 

1384 // — order the points in the same direction as XIO — > X20 
1385 

1386 if (XIO < X20) 

1387 

1388 { 
1389 

1390 xl = rootl ; 

1391 

1392 x2 = root2 ; 

1393 

1394 } 
1395 

1396 else 
1397 

1398 { 
1399 

1400 xl = root2 ; 

1401 

1402 x2 = rootl ; 

1403 

1404 } 

1405 

1406 

1407 

1408 // — The y— values can be calculated by substituting the 
1409 

1410 // — X— values into the line equation y — YIO -f- nvf(x — XIO) 
1411 

1412 yl = YIO + m*(xl - XIO) ; 

1413 

1414 y2 = YIO + m*(x2 - XIO) ; 

1415 

1416 } 
1417 

1418 else // — roots are y— values 

1419 

1420 { 
1421 

1422 // — order the points in the same direction as YIO — > Y20 
1423 

1424 if (YIO < Y20) 

1425 

1426 { 
1427 

1428 y 1 = root 1 ; 

1429 

1430 y2 = root2 ; 

1431 

1432 } 
1433 

1434 else 
1435 

1436 { 
1437 

1438 yl = root2 ; 

1439 

1440 y2 = rootl ; 

1441 

1442 } 

1443 

1444 

1445 

1446 // — The X— values can be calculated by substituting the 
1447 

1448 / / — y— V alues into the line equation x — XIO -f rt^ ( V — ) 
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1449 




1450 




1451 




1452 




1453 




1454 


} 


1455 




1456 




1457 




1458 


//- 


1459 




1460 


//- 


1461 




1462 


//- 


1463 




1464 


segi 


1465 




1466 




1467 




1468 


//- 


1469 




1470 


//- 


1471 




1472 


i f 


1473 




1474 


{ 


1475 




1476 




1477 




1478 


} 


1479 




1480 


C 1 S 1 


1481 




1482 


{ 


1483 




1484 




1485 




1486 




1487 




1488 


} 


1489 




1490 } 





xl ^ XIO + m*(yl - YIO) ; 
x2 ^ XIO + m*(y2 - YIO) ; 



McssagcCode ) : 



return — 1 ; 



(*Mcss age Code) ^ TWOJNTERSECTION_POINTS ; 
rcLurii segment_arca; 



Listing 14. C-SOURCE CODE FOR EL- 

LIPSE.ELLIPSE.OVERLAP 

6. APPENDIX C. 

1491 /* 

1492 
1493 * 
1494 

1495 ^- Function : double ellip s e ^ ellip s e ^ a v erl ap 

1496 

1497 * 

1498 

1499 * Purpose: Given the parameters of two ellipses . this function calculates 
1500 

1501 * the area of overlap between the two curves. If the ellipses 

are 

1502 

1503 * disjoint , this function returns 0.0; if one ellipse is 

contained 

1504 

1505 * within the other . this function returns the area of the 

enclosed 

1506 

1507 * ellipse ; if the ellipses intersect , this function returns the 

1508 
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1509 
1510 
1511 
1512 
1513 
1514 
1515 
1516 
1517 
1518 
1519 
1520 
1521 
1522 
1523 
1524 
1525 
1526 
1527 
1528 
1529 
1530 
1531 
1532 
1533 
1534 
1535 

1536 
1537 

1538 
1539 
1540 
1541 
1542 
1543 
1544 
1545 
1546 
1547 
1548 
1549 
1550 
1551 
1552 
1553 
1554 
1555 
1556 
1557 
1558 
1559 
1560 
1561 
1562 
1563 

1564 
1565 
1566 
1567 



calculated area of overlap. 

Reference: Hughes and Chraibi (2011), Calculating Ellipse Overlap Areas 

Dependencies : math.h for calls to trig and absolute value functions 
program^constants.h error message codes and constants 

Inputs : 



1 . 


double 


PHI_1 


GCW rotation angle 


of first ellipse , radians 


2. 


double 


Al 


semi— axis length in 


X— dire ction first ellipse 


3. 


double 


Bl 


semi— axis length in 


y— dire ction first ellipse 


4- 


double 


HI 


horizontal o ffs et o 


f center first ellipse 


5. 


double 


Kl 


vertical offset of 


center first ellipse 


6. 


double 


PHI.2 


CCW rotation angle 


of second ellipse , radians 


7. 


double 


A2 


semi— axis length in 


X— dire ction second 



ellipse 
ellipse 



Outputs : 



Return : 



8 . double B2 semi— axis length in y— dire ction second 

9 . double H2 horizontal a ffs e t of center second ellipse 
10. double K2 vertical offset of center second ellipse 

1 . int ^ rtnC o de returns diagnostic inforin ation integer code 
integ er codes in programme onstants . h 

The calculated value of the ov erlap area 

— 1 is returned in case of an error with the calculation 
is returned if the ellipses are disjoint 

pi^A^B of smaller ellipse if one ellipse is contained within 
the other ellipse 



// 



1568 

1569 //== DEFINE PROGRAM CONSTANTS 
1570 

1571 // 



1572 

1573 ^include " pro gr am_c o n s t an t s . h" // error message codes and constants 
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1574 
1575 
1576 

1577 // 



1578 

1579 //== DEPENDENT FUNCTIONS 



1580 
1581 



1582 

1583 double nointpts (douljlc Al , double Bl , double A2 , double B2 , double HI, 
1584 

1585 double Kl , double H2_TR , double K2_TR , double AA, double 

BB, 

1586 

1587 double CC, double DD, double EE, double FF , int *rtnCodo); 

1588 

1589 

1590 

1591 doul:)le twointpts (double xint [] , double yint [] , double Al , double Bl , 
1592 

1593 double PHI.l , double A2 , double B2 , double H2.TR , 

1594 

1595 double K2.TR , double PHI.2 , double AA, double BB, 

1596 

1597 double CC, double DD, double EE, double FF , int *rtnCodo) 

1598 
1599 
1600 

1601 double threeintpts (double xint [] , dotible yint [] , double Al , double Bl , 
1602 

1603 double PHI.l, double A2 , double B2 , double H2.TR , 

1604 

1605 double K2.TR , double PHI.2, double AA, double BB, 

1606 

1607 double CC, double DD, double EE, double FF , 

1608 

1609 int *rtnCode ) ; 

1610 

1611 

1612 

1613 double fourintpts (double xint [] , double yint [] , double Al , double Bl , 
1614 

1615 double PHI.l, double A2 , double B2 , double H2.TR , 

1616 

1617 double K2.TR , double PHI.2, double AA, double BB, 

1618 

1619 double CC, double DD, douljle EE, double FF , int *rtnCode 

) ; 

1620 
1621 
1622 

1623 int istanpt (double x, double y, double Al , dotible Bl , double AA, double BB 
1624 

1625 double CC, double DD, double EE, double FF) ; 

1626 

1627 

1628 

1629 douljle ellipsc2tr (double x, double y, double AA, double BB, 
1630 

1631 double CC, double DD, double EE, double FF) ; 

1632 

1633 

1634 

1635 // — functions for solving the quartic equation from Netlib/TOMS 
1636 

1637 void BIQUADROOTS (double p [] , double r[][5]); 



CALCULATING ELLIPSE OVERLAP AREAS 



49 



1638 

1639 void CUBICROOTS (double p [] , double r[][5]) 
1640 

1641 xoid QUADROOIS (double p[], double r[][5]); 

1642 

1643 

1644 

1645 // 



1646 
1647 



1648 
1649 



/== ELLIPSE-ELLIPSE OVERLAP 



1650 

1651 douljle e 1 1 i p s e _ c 1 1 i p s e _ o V e r 1 a p (double PHI_1 , double Al , double Bl , 
1652 

1653 double HI, double Kl , double PHI.2 , 

1654 
1655 
1656 
1657 
1658 
1659 { 
1660 
1661 



double A2, double B2 , double H2 , double K2 , 
int *rtnCode ) 



// 



1662 
1663 



1664 
1665 



//== DEFINE LOCAL VARIABLES 
// 



1666 
1667 
1668 
1669 
1670 
1671 
1672 
1673 
1674 
1675 
1676 
1677 
1678 
1679 
1680 
1681 
1682 
1683 
1684 
1685 



int i , j , k, nroots , nychk , nintpts , fnRtnCodc ; 

double AA, BB, CC, DD, EE, FF , H2.TR, K2.TR , A22 , B22 , PHI.2R; 

double cosphi , cosphi2 , sinphi , sinphi2 , cosphisinphi ; 

double tmpO , tmpl , tmp2 , tmp3 ; 

double cy[5] = {0.0}, py [ 5 ] = {0.0}, r[3][5] = {0.0}; 

double xl , x2 , yl2 , y22 ; 

double ychk[5] = {0.0}, xint[5], yint[5]; 

double Areal , Arca2 , OvcrlapArca; 



// 



1686 
1687 



1688 
1689 



//= DATA CHECK 

// 



1690 
1691 
1692 
1693 

1694 



// — Each of the ellipse axis lengths must he positive 

if ( (!(A1 > 0.0) \textbar \textbar ! ( Bl > 0.0)) \toxtbar \textbar (!( 
A2 > 0.0) \tcxtbar \tcxtbar ! ( B2 > 0.0)) ) 
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1695 


r 
I 


1696 




1697 


1 


1698 




1699 


1 


1700 




1701 


} 


1702 








1704 




1705 


// — 


1706 




1707 


if ( 


1708 




1709 


] 


1710 




1711 


if ( 


1712 




1713 


] 


1714 




1715 




1716 




1717 


// 


1718 




1719 


/ /== 


1720 




1721 


// 


1722 




1723 


//— 


1724 




1725 


// — 


1726 




1727 


// — 


1728 




1729 


// — 


1730 




1731 


// 


1732 




1733 




1734 




1735 


//— 


1736 




1737 


// — 


1738 




1739 


// — 


1740 




1741 


// — 


1742 




1743 


//— 


1744 




1745 




1746 




1747 




1748 




1749 


//— 


1750 




1751 


//— 


1752 




1753 


//— 


1754 




1755 


//— 


1756 




1757 


//— 



(*rtnCodc) ^ ERROR_ELLIPSE_PARAMETERS; 
return —1.0; 



PHI_1 ^ fmod {PHI_1, twopi): 



PHI_2 ^ fmod (PHI_2, twopi); 



//= DETEIUdlNE THE TWO ELLIPSE EQUATIONS FROM INPUT PARAMETERS 



Finding the p oints of intersection h etween two g eneral ellipses 

r e quir e s s olving a quartic equatioii . B efo re atteinpting to solve 
the 

quartic , several quick tests can he used to eliminate some cases 



at the problem , by addressing the easiest cases first. 

Working with the tr anslate d-j-r otated ellipses simplifies the 

calculations . The ellipses are translated then r otated so that 
the 

first ellipse is centered at the origin and oriented with the 
CO or din at e axes . Theii , the first ellipse will have the implicit 
( p olyno7nial ) form of 

x\^{}2/Al\^{}2 + y+2/Bl\^{}2 - 1 



Kl) 



/ / through the angle —PHI_1 . 
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1758 

1759 // The translated and rotated center point coordinates for the second 

1760 

1761 // ellipse are found with the rotation matrix, derivations are 

1762 

1763 // described in the reference. 

1764 

1765 cosphi = cos (PHI.l); 

1766 

1767 sinphi = sin (PHI.l); 

1768 

1769 H2.TR = (H2 - HI)* cosphi + (K2 - Kl)* sinphi; 

1770 

1771 K2.TR = (HI - H2)*sinphi + (K2 - Kl)*cosphi; 

1772 

1773 PHI.2R = PHI.2 - PHI.l ; 

1774 

1775 il ( (fabs (PHI.2R) > (twopi)) ) 

1776 

1777 PHI.2R = fmod (PHI.2R, twopi); 

1778 

1779 

1780 

1781 // Calculate implicit (Polynomial) coefficients for the second 

ellipse 

1782 

1783 // in its translated —by (—HI, —H2) and rotated— by —PHl^l postion 

1784 

1785 // — AA*x~{}2 + BB*x*y + CC*y~{}2 + DD* x + EE*y + FF = 

1786 

1787 // Formulas derived in the reference 

1788 

1789 // To speed things up, store multiply — used expressions first 

1790 

1791 cosphi = cos (PHI.2R) ; 

1792 

1793 cosphi2 — cosphi*cosphi; 

1794 

1795 sinphi = sin (PHI.2R) ; 

1796 

1797 sinphi2 — s i n p h i * s i n p h i ; 

1798 

1799 cosphisinphi — 2.0*cosphi*sinphi; 

1800 

1801 A22 = A2*A2; 

1802 

1803 B22 = B2*B2 ; 

1804 

1805 tmpO = ( cosphi *H2.TR + s i n p h i *K2.TR) / A22 ; 

1806 

1807 tnipl = ( sinphi *H2.TR - cosphi *K2.TR) /B22 ; 

1808 

1809 tnip2 = cosphi *H2.TR + s i n p h i *K2.TR ; 

1810 

1811 tnip3 = sinphi *H2.TR - c os p h i *K2.TR ; 

1812 

1813 

1814 

1815 // — implicit polynom.ia,l coefficients for the second ellipse 
1816 

1817 AA = cosphi2/A22 + sinphi2/B22; 

1818 

1819 BB = cosphisinphi/A22 — c o s p h i s i n p h i / B22 ; 

1820 

1821 CC = sinphi2/A22 + cosphi2/B22; 

1822 

1823 DD = —2.0* cosphi *tmpO — 2 . * s i n p h i * tmpl ; 

1824 

1825 EE = -2.0* sinphi *tmpO + 2 . * cosp h i * tmpl ; 

1826 

1827 FF = tmp2*tmp2/A22 + tmp3* tmp3/B22 - 1.0; 

1828 

1829 



52 



GARY B. HUGHES AND MOHCINE CHRAIBI 



1830 
1831 



// 



1832 
1833 



1834 
1835 



//== CREATE AND SOLVE THE QUARTIC EQUATION TO FIND INTERSECTION POINTS 



// 



1836 
1837 
1838 
1839 
1840 
1841 
1842 
1843 
1844 
1845 
1846 
1847 
1848 
1849 
1850 
1851 
1852 
1853 
1854 
1855 
1856 
1857 
1858 
1859 
1860 
1861 
1862 
1863 
1864 
1865 
1866 
1867 
1868 
1869 
1870 
1871 
1872 
1873 
1874 
1875 
1876 
1877 

1878 
1879 
1880 
1881 
1882 
1883 
1884 
1885 
1886 
1887 
1888 
1889 
1890 
1891 
1892 
1893 
1894 
1895 
1896 



// // execution arrives here . the ellipses are at least 'close ' to 

/ / intersecting . 

// Coefficients for the Quartic Polynomial in y are calculated from 

// the two implicit equations . 

// — Formulas for these coefficients are derived in the reference. 
cy[4] = pow (Al, 4.0)*AA*AA + B1*B1*(A1»A1»(BB*BB - 2.0*AA*CC) 

+ B1*B1*CC*CC) ; 
cy[3] = 2.0*B1*(B1*B1*CC*EE + A1*A1*(BB*DD - AA*EE) ) ; 
cy[2] = A1*A1*((B1*B1*(2.0*AA*CC - BB*BB) + DD*DD - 2.0*AA*FF) 

- 2.0*A1*A1*AA*AA) + B1*B1 * ( 2 . * CC*FF + EE*EE) ; 
cy[l] = 2.0*B1*(A1*A1*(AA*EE - BB*DD) + EE*FF) ; 
cy[0] = (A1*(A1*AA - DD) + FF) * (Al* (A1*AA + DD) + FF) ; 

// Once the coefficients for the Quartic Equation in y are known, the 

// roots of the quartic polynomial will represent y— values of the 

// intersection points of the two ellipse curves. 

// The quartic sometimes degenerates into a polynomial of lesser 

/ / degree . so handle all possible cases. 

if (fabs (cy[4]) > 0.0) 
{ 

//== QUARTIC COEFFICIENT NONZERO. USE QUARTIC FORMULA 



} 



for (i = 0; i <= 3; i++) 

Py[4- i ] = cy [ i ] / cy [4] : 
py[0] = 1.0; 

BIQUADROOTS ( py , r ) ; 
nroots — 4; 

;lsc if (fabs (cy[3]) > 0.0) 



{ 
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1897 




/ / DTI A TiTTr^ FMT'D T?ATT?T^ A TTTQ TD DTTT^T/^ TTQTT DTTPIT/^ JPf~t]^Iii/fT IT A 

/ / \^Ut\ix 1 ILy UH/ V iiiX/i 1 HiO 1 kJ U IJIK^ , Uotli O U Dl\^ r \Ji\lvL ULiri. 


1898 






1899 




lOl \\ — U, 1 <. — A . 1 i — h j 


1900 






1901 




r -1 r • 1 / r 1 
Pyl3- 1 J ^ cy [ 1 J / cy [ 3 J ; 


1902 






1903 




py[0] - 1.0; 


1904 






1905 






1906 






1907 






1908 






1909 




11 r o o t s — o J 


1910 






1911 


\ 
/ 




1912 






1913 


CSC 




1914 






1915 


r 
i 




1916 






1917 




/ / DTI A nTTr' DTT'D TPATTPT^ A TTT'Q TD DTI A DT? A TTD TIQTP DTI A OR A TTD T?DT^J<,/fT IT A 


1918 






1919 




lOi (_ 1 — U; 1 < — i; iH — V) 


1920 






1921 




Pyl2- 1 1 = cy [ 1 ] / cy [ 2 J ; 


1922 






1923 




py[0] = 1.0; 


1924 






1925 






1926 






1927 






1928 






1929 




nroots — z; 


1930 






1931 


\ 
} 




1932 






1933 


else 


11 (^laus i^cy[ijj > u.uj 


1934 






1935 


r 
1 




1936 






1937 




/ / DTI A TiTTD DTPDTT'ATTT'T? A TTPQ TD T TAIfTA P/ ■ QDT l/P* DT/? fPDTT V 


1938 






1939 




//— cy[l]* Y + cy [0] = 


1940 






1941 




r I -L J I -L J — I — cy i^J / cy [ i J J , 


1942 






1943 




1- r 9 1 r 1 1 n n ■ 


1944 






1945 




nroots — 1| 


1946 






1947 


\ 
/ 




1948 






1949 


else 




1950 






1951 


I 




1952 






1953 




//= COMPLETELY DEGENERATE QUARTIC: ELLIPSES IDENTICAL??? 


1954 






1955 




/ / — a completely deg enerate quartic , which would seem to 


1956 






1957 




/ / — indicate that the ellipses are identical. However . somt 


1958 






1959 




/ / — c o nfi gurations lead to a deg enerate guar tic with no 


1960 






1961 




// — points of inter sectio7i . 


1962 






1963 




nroots — 0; 


1964 






1965 


} 
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1966 
1967 
1968 

1969 // 
1970 

1971 //== CHECK ROOTS OF THE QUARTIC: ARE THEY POINTS OF INTERSECTION? 



1972 

1973 // 



1974 

1975 // determine which roots are real, discard any complex roots 

1976 

1977 nychk = 0; 

1978 

1979 lor (i = 1; i <= nroots; i++) 

1980 

1981 { 
1982 

1983 il ( fabs ( r [2] [ i ] ) < EPS) 

1984 

1985 { 
1986 

1987 nychk + + ; 

1988 

1989 ychk[ nychk] = r [ 1 ] [ i ] * Bl ; 

1990 

1991 } 
1992 

1993 } 

1994 

1995 

1996 

1997 // sort the real roots by straight insertion 

1998 

1999 for (j = 2; j <= rrychk ; j ++) 

2000 

2001 { 
2002 

2003 tmpO = ychk [ j ] ; 

2004 

2005 

2006 

2007 for (k = j- 1; k >= 1: k ) 

2008 

2009 { 
2010 

2011 if (ychk[k] <= tmpO) 

2012 

2013 break ; 

2014 

2015 

2016 

2017 ychk[k + l] = ychk[k]; 

2018 

2019 } 

2020 

2021 

2022 

2023 ychk[k + l] = tinpO ; 

2024 

2025 } 

2026 

2027 

2028 

2029 // determine whether polynomial roots are points of intersection 

2030 

2031 // — for the two ellipses 
2032 

2033 nintpts = 0; 
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2034 
2035 
2036 
2037 
2038 
2039 
2040 
2041 
2042 
2043 
2044 
2045 
2046 
2047 
2048 
2049 
2050 
2051 
2052 
2053 
2054 
2055 
2056 
2057 
2058 
2059 
2060 
2061 

2062 
2063 
2064 
2065 
2066 
2067 
2068 
2069 
2070 
2071 
2072 
2073 
2074 
2075 
2076 
2077 
2078 
2079 
2080 
2081 
2082 
2083 
2084 
2085 

2086 
2087 
2088 
2089 
2090 
2091 
2092 
2093 
2094 
2095 
2096 
2097 
2098 
2099 
2100 
2101 
2102 
2103 
2104 



for (i — 1; i <— nychk ; iH — h) 
{ 

// — check for multiple roots 

if ((i > 1) \&\& (fabs (ychk[i] - ychk[i-l]) < (EPS/2.0); 
continue ; 

// — check iriicrscclion points for ychkfi] 
if ( fabs (ychk[ i ] ) > Bl ) 
xl = 0.0; 

else 

xl = Al*sqrt (1.0 - (ychk [ i ] * ychk [ i ] ) /(B1*B1) ) ; 
x2 = -xl : 



if ( fabs ( ollipsc2tr (xl , ychk [ i ] , AA, BB , CC, DD, EE, FF) ) < EPS 
/2.0) 



n i n t p t s H — h; 

il' (niiitpts > 4) 



{ 



} 



(♦rtnCode) = ERRORJNTERSECTION JTS ; 
return — 1.0: 



xint[nintpts] — xl; 

y i n t [ n i n t p t s ] — ychk [ i 



if (( fabs ( ellipse2tr (x2 , ychk [ i ] , AA, BB, CC, DD. EE, FF) ) < EPS 
/2.0) 

\&\& (fabs (x2 - xl) > EPS/2.0)) 

r 

n i n t p t s H — h; 

if (nintpts > 4) 

{ 

(♦rtnCode) = ERRORJNTERSECTION JTS ; 
return —1.0; 



} 

xint [nintpts] — x2; 
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2105 y i n t [ n i n t p t s ] — ychk [ i ] ; 

2106 

2107 } 
2108 

2109 } 

2110 

2111 

2112 

2113 // 



2114 

2115 //== HANDLE ALL CASES FOR THE NUMBER OF INTERSCTION POINTS 



2116 

2117 // 



2118 

2119 switch (nintpts) 

2120 

2121 { 
2122 

2123 case 0: 

2124 

2125 case 1: 

2126 

2127 OverlapArea = nointpts (Al, Bl , A2 , B2 , HI, Kl , H2.TR , K2.TR , 

AA, 

2128 

2129 BB, CC, DD, EE, FF , rtnCodo); 

2130 

2131 return OverlapArea; 

2132 

2133 

2134 

2135 case 2: 

2136 

2137 // when there are two intersection points . it is possible for 

2138 

2139 // them to both be tangents , in which case one of the 

ellipses 

2140 

2141 // is fully contained within the other. Check the points for 

2142 

2143 // tang ents ; if one of the points is a tangent . then the 

other 

2144 

2145 // must be as well . otherwise there would be more than 2 

2146 

2147 // intersection points. 

2148 

2149 fnRtnCode = istanpt (xint[l], yint[l], Al , Bl , AA, BB, CC, DD, 

2150 

2151 EE, FF) : 

2152 

2153 

2154 

2155 if (fnRtnCode = TANGENT JOINT) 

2156 

2157 OverlapArea = nointpts (Al, Bl , A2 , B2 , HI, Kl , H2_TR, 

K2.TR , 

2158 

2159 AA, BB, CC, DD, EE, FF, rtnCode); 

2160 

2161 else 
2162 

2163 OverlapArea = twointpts (xint, yint, Al , Bl , PHI.l , A2 , B2 , 

2164 

2165 H2.TR, K2.TR, PHI.2 , AA, BB, CC, 

DD, 

2166 

2167 EE, FF, rtnCode) ; 
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2168 
2169 
2170 
2171 
2172 
2173 
2174 
2175 
2176 
2177 
2178 
2179 
2180 
2181 
2182 
2183 
2184 
2185 
2186 
2187 
2188 
2189 
2190 
2191 
2192 
2193 
2194 
2195 
2196 
2197 
2198 
2199 
2200 
2201 
2202 
2203 
2204 
2205 

2206 
2207 
2208 
2209 
2210 
2211 
2212 
2213 
2214 
2215 
2216 
2217 
2218 
2219 
2220 
2221 
2222 
2223 
2224 
2225 
2226 
2227 
2228 
2229 
2230 
2231 
2232 
2233 

2234 
2235 
2236 
2237 
2238 



return OvcrlapArca; 
case 3 : 

// when there are three intersection points , one and only one 

// of the points must be a tangent point. 

OverlapArea = threeintpts (xint, yint , Al , Bl , PHI.l , A2 , B2 , 

H2.TR, K2.TR, PHI.2 , AA, BE, CC, DD, 
EE, FF, rtnCode) ; 

return OverlapArea; 
case 4 : 

// four intersections points has only one case. 

OverlapArea = fourintpts (xint, yint, Al , Bl , PHI.l, A2 , B2 , 

H2.TR, K2.TR, PHI.2, AA, BB, CC, DD, 

EE, FF, rtnCode ) ; 

return OverlapArea; 
default : 

// should never get here (but get compiler warning for 

missing 

// return value if this line is omitted) 

(*rtnCodc) = ERRORJNTERSECTION JTS ; 
return — 1.0; 



double ellipse2tr (double x, double y, double AA, double BB, 

double CC, double DD, double EE, double FF) 

return (AA*x*x + BB*x*y + CC*y*y + DD*x + EE*y + FF) ; 



{ 



doul-jlc nointpts (double Al , double Bl , double A2 , double B2 , double HI, 

double Kl, double H2.TR , double K2.TR , double AA, double 
BB, 

double CC, double DD, double EE, double FF , int *rtnCode) 
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2239 
2240 
2241 
2242 
2243 
2244 
2245 
2246 
2247 
2248 
2249 
2250 
2251 
2252 
2253 
2254 
2255 
2256 
2257 
2258 
2259 
2260 
2261 
2262 
2263 
2264 
2265 
2266 
2267 
2268 
2269 
2270 
2271 
2272 
2273 
2274 
2275 
2276 
2277 
2278 
2279 
2280 
2281 
2282 
2283 
2284 
2285 
2286 
2287 
2288 
2289 
2290 
2291 
2292 
2293 
2294 
2295 
2296 
2297 
2298 
2299 
2300 
2301 
2302 
2303 
2304 
2305 
2306 
2307 
2308 
2309 
2310 
2311 



// The relative size of the two ellipses can be found from the axis 

/ / lengths 

double relsizo = (A1*B1) - (A2*B2); 

if (rolsizc > 0.0) 
{ 

/ / — First Ellipse is larger than second ellipse. 
// — // second ellipse center (H2^TR. K2^TR) is inside 
/ / — first ellipse , then ellipse 2 is completely inside 
/ / — ellipse 1 . Otherwise . the ellipses are disjoint. 
if ( ((H2.TR*H2.TR) / (A1*A1) 

+ (K2.TR*K2.TR) / (B1*B1)) < 1.0 ) 

{ 

(*rtnCodc) = ELLIPSE2.INSIDE.ELLIPSE1 ; 
return ( pi *A2*B2) ; 



} 



J 1 s c 



(*rtnCodc) ^ DISJOINT.ELLIPSES ; 
r c t Li r 11 0.0: 



} 



:lsc if (rclsizc < 0.0) 



{ 



/ / — Second Ellipse is larger than first ellipse 

/ / — // first ellipse center ( . 0) is inside the 

/ / — second ellipse, then ellipse 1 is completely inside 

/ / — ellipse 2 . Otherwise , the ellipses are disjoint 

// — A4*x"{}S + BB^x^y + CC^y\^{}2 + DD^x -h EE^y + FF ^ 

if (FF < 0.0) 



{ 



(*rtnCodc) ^ ELLIPSE1_INSIDE_ELLIPSE2 ; 
return (pi*Al*Bl) ; 



? 1 s c 



(*rtnCodc) ^ DISJOINT.ELLIPSES ; 
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2312 

2313 return 0.0; 

2314 

2315 } 
2316 

2317 } 
2318 

2319 else 
2320 

2321 { 
2322 

2323 // — // execution arrives here, the relative sizes are identical. 
2324 

2325 // — Are the ellipses the sa/nie? Check the parameters to see. 
2326 

2327 if ((HI = H2.TR) \&\& (Kl = K2.TR) ) 

2328 

2329 { 
2330 

2331 (*rtnCode) = ELLIPSES_ARE JDENTICAL : 

2332 

2333 return (pi*Al*Bl); 

2334 

2335 } 
2336 

2337 else 
2338 

2339 { 
2340 

2341 // should never get here . so return error 

2342 

2343 (*rtnCode) = ERROR-CALCULATIONS; 

2344 

2345 return -1.0; 

2346 

2347 } 
2348 

2349 } // end if (r el size > 0.0) 

2350 
2351 } 
2352 
2353 
2354 

2355 // — two distinct intersection points (xl, yl) and (x2 , y2) find overlap 
area 

2356 

2357 doul:)lc twointpts (double x[] , double y[] , double Al , double Bl , double 
PHI.l , 

2358 

2359 double A2 , double B2 , double H2.TR , double K2.TR , 

2360 

2361 double PHI.2 , double AA, double BB, double CC, double DD, 

2362 

2363 double EE, double FF , int *rtnCode) 

2364 

2365 { 

2366 

2367 double areal , area2 ; 

2368 

2369 double xmid , ymid , xmid_rt , ymid_rt ; 

2370 

2371 double thetal, thcta2; 

2372 

2373 double tmp , trsign ; 

2374 

2375 double xl.tr, yl.tr, x2_tr, y2_tr; 

2376 

2377 double discr; 

2378 

2379 double cosphi , sinphi ; 

2380 

2381 

2382 
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2383 // if execution arrives here, the intersection points are not 

2384 

2385 // tangents . 

2386 
2387 
2388 

2389 // determine which direction to integrate in the ellipse-segment 

2390 

2391 // routine for each ellipse. 

2392 
2393 
2394 

2395 // find the parametric angles for each point on ellipse 1 

2396 

2397 if (fabs (x[l]) > Al) 

2398 

2399 x[l] = (x[l] < 0) ? -Al : Al ; 

2400 

2401 if (y[l] < 0.0) // — QuadranL III or IV 

2402 

2403 thetal = twopi - acos (x[l] / Al ) ; 

2404 

2405 else /'/ — Quadrant I or II 

2406 

2407 thetal = acos (x[l] / Al); 

2408 

2409 

2410 

2411 if (fabs (x[2]) > Al ) 

2412 

2413 x[2] = (x[2] < 0) ? -Al : Al ; 

2414 

2415 if (y[2] < 0.0) // — Quadrant III or IV 

2416 

2417 thGta2 = twopi — acos (x[2] / Al ) ; 

2418 

2419 else // — Quadrant I or II 

2420 

2421 theta2 = acos (x[2] / Al); 

2422 

2423 

2424 

2425 // — logic is for proceeding counterclockwise from thetal to theta2 
2426 

2427 if (thetal > thota2) 

2428 

2429 { 
2430 

2431 tmp = thetal ; 

2432 

2433 thetal = thcta2; 

2434 

2435 theta2 = tmp; 

2436 

2437 } 

2438 

2439 

2440 

2441 // find a point on the first ellipse that is different than the two 

2442 

2443 // intersection points. 

2444 

2445 xmid = Al*cos ((thetal + theta2)/2.0); 

2446 

2447 ymid = Bl*sin ((thetal + theta2)/2.0); 

2448 
2449 
2450 

2451 // the point (xmid, ymid) is on the first ellipse ^between' the two 

2452 

2453 // intersection points (xflj, yflj) and (x[2] , y[2]) when travelling 

2454 
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2455 // counter— clockwise from ( x [ 1 ] . y [ 1 ] ) to (x [ 2] , y [ 2] ) . If the 

point 

2456 

2457 // (xmid, ymid) is inside the second ellipse , then the desired 

segment 

2458 

2459 // of ellipse 1 contains the point (xmid, ymid) , so integrate 

2460 

2461 / / counterclockwise from (x [ 1 ] . y [ 1 ] ) to ( x [ 2] , y [ 2] ) . Otherwise , 

2462 

2463 // integrate counterclockwise from (x[2] . y[2j) to (xflj, y[l]) 

2464 

2465 if (cllipso2tr (xmid, ymid, AA, BB, CC, DD, EE, FF) > 0.0) 

2466 

2467 { 
2468 

2469 tmp = thetal ; 

2470 

2471 thetal = thcta2 ; 

2472 

2473 thcta2 = tmp; 

2474 

2475 } 

2476 

2477 

2478 

2479 // here is the ellipse segment routine for the first ellipse 

2480 

2481 if (thetal > theta2) 

2482 

2483 thetal — = twopi; 

2484 

2485 if ((theta2 - thetal) > pi) 

2486 

2487 trsign = 1.0; 

2488 

2489 else 
2490 

2491 trsign = -1.0; 

2492 

2493 areal = . 5 * ( A1*B1 * ( t ho t a 2 - thetal) 

2494 

2495 + trsign*fabs (x[l]*y[2] - x[2]*y[l])); 

2496 

2497 

2498 

2499 // find ellipse 2 segment area. The ellipse segment routine 

2500 

2501 // needs an ellipse that is centered at the origin and oriented 

2502 

2503 // with the coordinate axes. The intersection points (x[l] . yfl]) 

and 

2504 

2505 // (x[2j, y[2j) are found with both ellipses translated and rotated 

by 

2506 

2507 // (—HI, —Kl) and —PHI^l. Further translate and rotate the points 

2508 

2509 // to put the second ellipse at the origin and oriented with the 

2510 

2511 // coordinate axes. The translation is (—H2-TR, —K2-TR), and the 

2512 

2513 // rotation is -(PHI.2 - PHI.l) = PHI.l - PHI.2 

2514 

2515 cosphi = cos (PHI.l - PHI. 2); 

2516 

2517 sinphi = sin (PHI.l - PHI. 2); 

2518 

2519 xl.tr = (x[l] - H2.TR) * eosphi + (y[l] - K2.TR)*-sinphi ; 

2520 

2521 yl.tr = (x[l] - H2.TR) * sinphi + (y[l] - K2.TR) * cosphi ; 

2522 

2523 x2.tr = (x[2] - H2.TR) * eosphi + (y[2] - K2.TR) *- s i n p h i ; 
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2524 
2525 
2526 
2527 
2528 
2529 
2530 
2531 
2532 
2533 
2534 
2535 
2536 
2537 
2538 
2539 
2540 
2541 
2542 
2543 
2544 
2545 
2546 
2547 
2548 
2549 
2550 
2551 
2552 
2553 
2554 
2555 
2556 
2557 
2558 
2559 
2560 
2561 
2562 
2563 
2564 
2565 
2566 
2567 
2568 
2569 
2570 
2571 
2572 
2573 
2574 
2575 
2576 
2577 
2578 
2579 
2580 
2581 
2582 
2583 
2584 
2585 
2586 
2587 
2588 
2589 
2590 
2591 
2592 
2593 
2594 
2595 
2596 



y2.tr = (x[2] - H2.TR) * s i n p h i + (y[2] - K2.TR) * cosphi ; 

// determine which branch of the ellipse to integrate by finding a 

// point on the second ellipse , and asking whether it is inside the 

// first ellipse ( in their once— tr an s I at e d+r ot at e d positions ) 

// — find the parametric angles for each point on ellipse 1 
if (fabs (xl.tr) > A2) 

xl.tr = (xl.tr < 0) ? -A2 : A2 ; 
if (yl.tr < 0.0) // — Quadrant III or IV 

thctal — twopi — acos (xl.tr/A2) ; 
else // — Quadrant I or II 

thctal — acos (xl.tr/A2); 

if (fabs (x2.tr) > A2) 

x2.tr = (x2.tr < 0) ? -A2 : A2 ; 
if (y2.tr < 0.0) // — Quadrant III or IV 

theta2 — twopi — acos (x2.tr/A2) ; 
else // — Quadrant I or II 

thcta2 = acos (x2.tr/A2); 

// logic is for proceeding counterclockwise from thctal to theta2 

if (thctal > theta2) 
{ 

tmp — thctal; 
thctal — thcta2; 
t hct a2 — tmp ; 

} 

/ / find a p oint on the second ellipse that is different than the two 

/ / intersection p oints . 

xmid ^ A2*cos ((thetal + theta2)/2.0); 
ymid ^ B2*sin ((thctal + t hct a2 ) / 2 . ) ; 

/ / translate the p oint back to the second ellipse in its one e— 

/ / translate d-hro tated position 

cosphi ^ cos (PHI_2 - PHI_1); 
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2597 
2598 
2599 
2600 
2601 
2602 
2603 
2604 
2605 

2606 
2607 
2608 
2609 
2610 
2611 
2612 
2613 
2614 
2615 
2616 
2617 
2618 
2619 
2620 
2621 
2622 
2623 
2624 
2625 
2626 
2627 
2628 
2629 
2630 
2631 
2632 
2633 
2634 
2635 
2636 
2637 
2638 
2639 
2640 
2641 
2642 
2643 
2644 
2645 
2646 
2647 
2648 
2649 
2650 
2651 
2652 
2653 
2654 
2655 
2656 
2657 
2658 
2659 
2660 
2661 
2662 
2663 
2664 
2665 
2666 
2667 
2668 



sinphi = sin (PHI. 2 - PHI.l); 

xmid_rt — xinid*cosphi + y mid*— s i n p h i + H2_TR ; 
ymid_rt — xmid* s i n p li i + ymid* c o s p li i + K2_TR ; 

// the point (xmid^rt, ymid^rt) is on the second ellipse 'between' 

the 

// intersection points (xflj, y[l]) and (x[2j . y[2j) when travelling 

/ / counterclockwise from (x [ 1 ] . y [ 1 j } to ( x [ 2] , y [ 2] ) . If the p oint 

// (xm,id_rt. ymid_rt) is inside the first ellipse , then the desired 

/ / segment of ellipse 2 contains the p oint ( xmid^rt , ymid^rt ) , so 

/ / integrate counterclockwise from ( x [ 1 J , y [ i ] ) to ( x [2 j . y [ 2] ) . 

/ / Otherwise , integrate counterclockwise from ( x [ 2] . y [ 2] ) to 

//— (x[l], y[l]) 

if (((xmid_rt*xmid_rt) /(A1*A1) + ( y mid.rt * y mid.rt ) / ( Bl * Bl ) ) > 1.0) 
{ 

tmp — t hot a 1 ; 
thctal — thcta2; 
t het a2 — tmp ; 



/ / here is the ellipse segment routine for the second ellipse 

if (thctal > thcta2) 

thctal — — two pi ; 
il' ({thcta2 - thctal) > pi) 

trsign — 1.0; 

else 

trsign — — 1.0; 
arca2 ^ . 5 * ( A2*B2 * ( t he t a 2 - thetal) 

+ trsign^fabs (xl_tr*y2_tr — x2_tr*yl_tr)); 

(*rtnCode) ^ TWOJNTERSECTION_POINTS ; 
return a r c a 1 + a r c a 2 ; 



/ / — three distinct intersection points , must have two inters actions 
/ / — and one tangent . which is the only possibility 

double thrccintpts (double xint [] , double yint [] , double Al , double Bl , 
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2669 

2670 

2671 

2672 

2673 

2674 

2675 

2676 

2677 { 

2678 

2679 

2680 

2681 

2682 

2683 

2684 

2685 

2686 

2687 

2688 

2689 

2690 

2691 

2692 

2693 

2694 

2695 

2696 

2697 

2698 

2699 

2700 

2701 

2702 

2703 

2704 

2705 

2706 

2707 

2708 

2709 

2710 

2711 

2712 

2713 

2714 

2715 

2716 

2717 

2718 

2719 

2720 

2721 

2722 

2723 

2724 

2725 

2726 

2727 

2728 

2729 

2730 

2731 

2732 

2733 

2734 

2735 

2736 

2737 

2738 

2739 

2740 

2741 



double PHI.l , double A2 , double B2 , double H2.TR , 
double K2.TR, double PHI.2 , double AA, double BB, 
double CC, double DD, double EE, double FF , 
int *rtnCQde) 

int i , tanpts , tanindcx , fnRtn ; 
double OverlapArea; 

// need to determine which point is a tangent , and which two points 

/ / — arc intersections 
tanpts — 0; 

for (i = 1; i <= 3; i++) 
{ 

fnRtn = istanpt (xint[i], yint[i], Al , Bl , AA, BB, CC, DD, EE, FF) : 



{ 



if (fnRtn = TANGENTJ^OINT) 

tanpts + + ; 
tanindcx — i: 



// there MUST he 2 intersection points and only one tang ent 

ii (tanpts !— 1) 



{ 



/ /— — should never get here unless there is a prohlein dis c erning 
/ / — whether or not a point is a tang ent or intersection 
(*rtnCode) ^ ERRORJNTERSECTION_PTS ; 
r c L n r 11 —1.0; 



/ / store the two inter es ection points into (x [ 1 ] , y [ 1 ] ) and 

//— (x[2], y[2]) 
switch (tanindcx) 

{ 

case 1 : 

xint[l] — xint[3]; 
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2742 
2743 
2744 
2745 
2746 
2747 
2748 
2749 
2750 
2751 
2752 
2753 
2754 
2755 
2756 
2757 
2758 
2759 
2760 
2761 
2762 
2763 
2764 
2765 
2766 
2767 
2768 
2769 

2770 
2771 
2772 
2773 
2774 
2775 
2776 
2777 
2778 
2779 
2780 
2781 
2782 
2783 
2784 
2785 
2786 
2787 
2788 
2789 

2790 
2791 
2792 
2793 
2794 
2795 
2796 
2797 
2798 
2799 
2800 
2801 
2802 
2803 
2804 
2805 
2806 
2807 
2808 
2809 
2810 
2811 
2812 



yint [1] = yint [3] ; 
break ; 



case 2 : 

xint [2] 
yint [2] 
break ; 



— xint [ 3 ] 
= yint [3] 



case 3 : 

/ / intersection points are already in the right places 

break ; 



OverlapArea = twointpts (xint, yint, Al , Bl , PHI.l , A2 , B2 , H2_TR, 
K2.TR , 

PHI.2 , AA, BE, CC, DD, EE, FF , rtnCode ) : 
(*rtnCodo) = THREE JNTERSECTION J'OINTS ; 
return OverlapArea; 



/ / — four intersection p oints 

double fourintpts (double xint [] , double yint [] , double Al , double Bl , 
double PHI.l, double A2 , double B2 , double H2.TR , 
double K2.TR, double PHI.2 , double AA, double BB, 



double CC, double DD, double EE, double FF , int *rtnCode 
) 



int. i , j , k ; 

double xmid , ymid , xint_tr [5] , yint_tr [5] , OverlapArea; 

double theta[5] , theta_tr [5] , cosphi , sinphi , tmpO , tmpl , tmp2 ; 

double areal , area2 , areaS , area4 , arcaS ; 

//— only one case , which involves two segments from each ellipse , plus 
//— two triangles . 

/ / get the parametric angles along the first ellipse for each of the 

/ / intersection p oints 

for ( i ^ 1; i <^ 4; i++) 
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2813 
2814 
2815 
2816 
2817 
2818 
2819 
2820 
2821 
2822 
2823 
2824 
2825 
2826 
2827 
2828 
2829 
2830 
2831 
2832 
2833 
2834 
2835 
2836 
2837 
2838 
2839 
2840 
2841 
2842 
2843 
2844 
2845 
2846 
2847 
2848 
2849 
2850 
2851 
2852 
2853 
2854 
2855 
2856 
2857 
2858 
2859 
2860 
2861 
2862 
2863 
2864 
2865 
2866 
2867 
2868 
2869 
2870 
2871 
2872 
2873 
2874 
2875 
2876 
2877 
2878 
2879 
2880 
2881 
2882 
2883 
2884 
2885 



i r ( fabs ( xint [ i ] ) > Al) 

xint[i] = (xint[i] < 0) ? -Al : Al ; 
if (yint[i] < 0.0) // — Quadrant III or IV 

thcta[i] — twopi — acos (xint[i] / Al) ; 
else // Quadrant I or II 

thcta[i] — acos {xint[i] / Al); 



// sort the angles by straight insertion , and put the points in 

// e ounter— e I o ckiuis e order 

for (j = 2; j <= 4; j ++) 
{ 

tmpO — t liet a [ j ] ; 
tmpl — xint [ j ] ; 
tmp2 = y int [ j ] ; 



for (k = j - 1; k >= 1: k- 



{ 



if (theta [k] <= tmpO) 
break ; 

tiicta [k+1] = tlicta [k] ; 
xint [k + 1] = xint [k] ; 
yint [k + 1] = yint [k] ; 



tlieta [k+1] = tmpO; 
xint [k+1] — tmpl; 
yint [k+1] = tmp2 ; 



// find the area of the interior quadrilateral 

arcal = 0.5*fabs ((xint[3] - x i n t [ 1 ] ) * ( y i n t [ 4 ] - yint[2]) 

- (xint[4] - xint [2] ) *(yint [3] - yint[l])); 

// the intersection points lie on the second ellipse in its once 
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2887 


//— 


tr a7is late d-hrot ate d position . The segment alg orithin is iinplemented 


2888 
2889 


//— 


for an ellipse that is centered at the origin , and oriented with 


2890 






2891 


//— 


the coordiiiate axes ; so , in order to use the segment alg orithm 


2892 








//— 


with the second ellipse , the intersection p oiiits must he furth e r 


2894 






2895 


//— 


tr ans late d-hrot ate d by a'mounts that put the second ellipse c enter ed 


2896 






2897 


//— 


at the origin and oriented with the coordinate axes . 


2898 






2899 


cosphi ^ cos (PHI_1 - PHI_2); 


2900 






2901 


sinphi ^ sin (PHI_1 - PHI_2); 


2902 






2903 


for 


(i ^ 1; i <^ 4; i++) 


2904 






2905 


{ 




2906 






2907 




xint_tr[i] — (xint[i] — H2_TR) * cosphi + (yiiitfi] — K2_TR) *— s i n p h i ; 


2908 






2909 




yint_tr[ij — (^xint[ij — llz_l llj*sinpni + (^yiiit[ij — lv^_l rv. j * cospni ; 


2910 






291 1 






2912 






2913 




11 (tabs (xint_tr[iJJ > A2 ) 


2914 






2915 




xint_tr[i] — (xint_tr[i] < 0) ? — A2 : A2 ; 


2916 






2917 




if (yint_tr[i] < 0.0) // — Quadrant III or IV 


2918 






2919 




theta_ti'[ij — twopi — acos (xint_tr[ij / A2 ) ; 


2920 






2921 




else / /- — - Quadrant I or II 


2922 






2923 




thcta_tr[i] — acos (xint_tr[i] / A2 ) ; 


2924 






2925 


} 




2926 






2927 












2929 


//— 


g et the area of the two segments on ellipse 1 


2930 






2931 


xmid 


^ Al*cos ((thetafl] + t het a [ 2 ] ) / 2 . ) ; 


2932 






2933 


ymid 


^Bl*sin ((thctafl] +theta[2])/2.0); 


2934 






2935 






2936 






2937 


//— 


the p oint (xmid , ymid) is on the first ellipse 'between ' the two 


2938 






2939 


//— 


sorted intersection p oints ( xint [ 1 J , yintfl]) and ( xint [ 2] , yint 






[ZD 


2940 






2941 


//— 


when travelling counter— do ckwis e from ( xint [ 1 ] , yintfl]) to 


2942 






2943 


//— 


( xint [ 2 ] , yint[2]). If the p oint ( xmid , ymid ) is inside the 






second 


2944 






2945 


//— 


ellipse , then one desired segment of ellipse 1 contains the point 


2946 






2947 


//— 


( xmid , ymid ) , so integrate counter do ckwis e from ( xint [ 1 ] , yint 






[in 


2948 






2949 


//— 


to ( xint [ 2 j . yint [2] ) for the first segment . and fro'm 


2950 






2951 


//— 


( xint [ 3 ] , yint[3] to ( xint [Jf. ] . yint [4 ] ) for the second segment . 


2952 






2953 


if ( 


Gllipse2tr (xmid, ymid, AA, BB, CC, DD, EE, EE) < 0.0) 


2954 






2955 


{ 
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2956 

2957 

2958 

2959 

2960 

2961 

2962 

2963 

2964 

2965 

2966 

2967 

2968 

2969 

2970 

2971 

2972 

2973 

2974 

2975 

2976 

2977 

2978 

2979 

2980 

2981 

2982 

2983 

2984 

2985 

2986 

2987 

2988 

2989 

2990 

2991 

2992 

2993 

2994 

2995 

2996 

2997 

2998 

2999 

3000 

3001 

3002 

3003 

3004 

3005 } 

3006 

3007 

3008 

3009 //- 
3010 

3011 int 

3012 

3013 

3014 

3015 { 

3016 

3017 

3018 

3019 

3020 

3021 

3022 

3023 

3024 
3025 
3026 



aroa2 = . 5 * ( Al *B1 * ( t he t a [ 2 ] - thcta[l]) 

- fabs ( xint [1] * yint [2] - xi n t [ 2 ] * y i n t [ 1 ] ) ) ; 
aroa3 = . 5 * ( Al *B1 * ( t he t a [ 4 ] - thcta[3]) 

- fabs ( xint [3] * yint [4] - x i n t [ 4 ] * y i n t [ 3 ] ) ) ; 
area4 = . 5 * ( A2*B2 * ( t h c t a _ t r [ 3 ] - thGta_tr[2]) 

— fabs (xint_tr[2]*yint_tr[3] — xint_tr[3]*yint_tr 
areas = . 5 * ( A2*B2 * ( t h e t a _ t r [ 1 ] - (theta_tr[4] - twopi)) 

— fabs (xint_tr[4]*yint_tr[l] — xint_tr[l]*yint_tr 



2]) 



4])) 



area2 = . 5 * ( Al *B1 * ( t he t a [ 3 ] - thcta[2]) 

- fabs ( xint [2] * yint [3] - x i n t [ 3 ] * y i n t [ 2 ] ) ) ; 
area3 = . 5 * ( Al *B1 * ( t he t a [ 1 ] - (thcta[4] - twopi)) 

- fabs (xint [4] * yint [1] - xi n t [ 1 ] * y i n t [ 4 ] ) ) ; 
area4 = . 5 * ( A2*B2 * ( t het a [ 2 ] - thcta[l]) 

— fabs (xint_tr[l]*yint_tr[2] — xint_tr[2]*yint_tr 
areas = . 5 * ( A2*B2 * ( t het a [ 4 ] - thcta[3]) 

— fabs (xint_tr[3]*yint_tr[4] — xint_tr[4]*yint_tr 



1])) 



[3]) 



OverlapArea — areal + area2 + area3 + area4 + areaS; 
(*rtnCodc) ^ FOURJNTERSECTION_POINTS ; 
1 c 1 1 1 1' M OverlapArea; 



check whether an intersection p oint is a tangent or a cross — point 
istanpt (double x, double y, double Al, double Bl, double AA, double BB 

double CC, double DD, double EE, double FF) 
double xl, yl, x2, y2, theta, testl , test2 , branch, eps_radian; 

/ / Avoid inv erse trig calculation errors : there could be an error 

/ / if \ t extb ar xl /A\ textbar > 1.0 wheii c ailing acos ( ) . If execution 

arriv e s here . 

/ / then the p oint is on the ellipse within EPS . 
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3027 


11 ^^laDs ) ^ ■'^^ ) 


3028 




3029 


X — (^x<U} : — Al : Ai ; 


3030 




3031 




3032 




3033 


/ / Cj alculate the. p cir cim etvic ang Le on the ellipse jot ( x , y J 


3034 




3035 


// 1 he p cir ciTfietric angles depend on the QUddrant where each point 






3037 


/ /— IS located. o ee 1 ab Le 1 in the r e j er en c e . 


3038 




3039 


ji ^^y<.u.uj / / \^ uaarant 1 1 1 ov i v 


3040 




3041 


tncta — twopi — acos (x / Al) ; 


3042 




3043 


else / /— Quadra'rit 1 or II 


3044 




3045 


tneta — acos (x / Alj ; 


3046 




3047 




3048 




3049 


/ / d e t e rriiin e the d. i s t a n c e jrom the origin to the point ( x , y ) 


3050 




3051 


brancii — sqrt (x*x -|- y^y) | 


3052 




3053 




3054 




3055 


/ /— — use the dist an c e to jmd a small angle , such that the distance 


3056 




3057 


// along ellipse 1 is approximately /i^ hiro 


3058 




3059 


;r /ViT-'inr.Vi inn n-j. t?T3C \ 
ji ^^Drancii <. ±uu.u* H/r^o ) 


3060 




3061 


cps_iaaian — z.u* rLr^o ; 


3062 




3063 


else 


3064 




3065 


eps_radian — asin ( 2 . * h^Po/ brancn ) ; 


3066 




3067 




3068 




3069 


/ /— — determine two p oints th at are on each side oj ( x , y ) and Lie on 


3070 




3071 


// the first ellipse 






3073 


xi — Al*cos (^tlicta --|~ eps_raaiaiij ; 


3074 




3075 


yl — Jjl^siii (^tiicta -f- eps_raaianj ; 


3076 




3077 


x2 — Al^cos (tncta — eps_radian) ; 


3078 




3079 


y2 — J31*siii (tncta — eps_radian); 


3080 




3081 




3082 




3083 


// — evaluate the two adjacent points in the second ellipse equation 






3085 


tcstl ^ ellipse2tr (xl , yl , AA, BB, CC, DD, EE, FF) ; 


3086 




3087 


test2 ^ cllipsc2tr (x2 , y2 , AA, BB, CC, DD, EE, FF) ; 


3088 




3089 




3090 




3091 


/ / if the ellipses are tang ent at the intersection point , then 


3092 




3093 


/ / p oints on h oth sides will either h oth he inside ellipse 1 , or 


3094 




3095 


/ / — they will h oth be outside ellipse 1 


3096 




3097 


ii ((testl*test2) > 0.0) 


3098 




3099 


r c t u 1 n TANGENT_POINT ; 
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3100 

3101 

3102 

3103 

3104 

3105 } 

3106 

3107 

3108 

3109 // 



else 



return INTERSECTION.POINT ; 



3110 

3111 //- 
3112 

3113 //- 
3114 

3115 //- 
3116 

3117 //- 
3118 

3119 //- 
3120 

3121 //- 
3122 

3123 //- 
3124 

3125 //- 
3126 

3127 //- 
3128 

3129 // 



CACM Algorithm 326: Roots of low order polynomials. 

Nonweiler , Terence R.F. , CACM Algorithm 326: Roots of low order 

polynomials , Communications of the ACM, vol. 11 no. 4, pages 

269-270 (1068). Translated into c and programmed by M. Dow. ANUSF, 

Australian National University , Canberra. Australia. 

Accessed at http : / / wurw . netlib . org/toms/326. 

Modified to void functions . integers replaced with floating p oint 
where appropriate , some other slight mo difications for readability 
and debugging ease. 



3130 
3131 
3132 
3133 
3134 
3135 
3136 
3137 
3138 
3139 
3140 
3141 
3142 
3143 
3144 
3145 
3146 
3147 
3148 
3149 
3150 
3151 
3152 
3153 
3154 
3155 
3156 
3157 
3158 
3159 
3160 
3161 
3162 
3163 
3164 
3165 
3166 
3167 
3168 



void QUADROOTS (double p[], double r[][5]) 



{ 



A 

Array r[3][5] p[5] 

Roots of poly p[0]* x'{}2 + p[l]*x + p[2] = 

x=r[l][k] + i r[2][k] k = l,2 

*/ 

double b , c , d ; 

b=-p[l]/(2.0*p[0]) ; 

c=p[2]/p[0]; 

d=b*b-c ; 

i f (d> = 0.0) 

{ 

il' (b>0.0) 

b=(r[l][2] = (sqrt(d)+b)); 

else 

b=(r[l][2] = (-sqrt(d)+b)) ; 
r[l][l] = c/b; 
r[2][l] = (r[2][2]=0.0) ; 
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3169 

3170 

3171 

3172 

3173 

3174 

3175 

3176 

3177 

3178 

3179 

3180 

3181 

3182 

3183 

3184 

3185 } 

3186 

3187 

3188 

3189 VI. 

3190 

3191 { 

3192 

3193 

3194 

3195 

3196 

3197 

3198 

3199 

3200 

3201 

3202 

3203 

3204 

3205 

3206 

3207 

3208 

3209 

3210 

3211 

3212 

3213 

3214 

3215 

3216 

3217 

3218 

3219 

3220 

3221 

3222 

3223 

3224 

3225 

3226 

3227 

3228 

3229 

3230 

3231 

3232 

3233 

3234 

3235 

3236 

3237 

3238 

3239 

3240 

3241 



} 

else 

{ 



d = (r[2][l]=sqrt(-d)); 

r[2][2]=-d; 

r[l][l] = (r[l][2]=b); 



id CUBICROOTS( double p[], double r[][5]) 

A 

Array r[3][5] p[5] 

R,oots of poly p [Oj* x\-{}3 + p[l]* x\-{}2 + p[2]*x + p[3j = 
x=r[l][kj + i r[2][k] k=l,...,3 
Assumes 0< ar dan ( x ) < p i /2 for x>0 

*/ 

double s,t,b,c,d; 
int. k ; 

If (p[0]! = 1.0) 
{ 

for (k=l;k<4;k++) 

p[k]=p[k]/p[0]; 
p[0]=1.0; 

} 

s=p [ 1 ] / 3 . ; 
t— s * p [ 1 ] ; 

b = 0.5*(s *( t/1.5-p [2] )+p [3] ) ; 
t=(t-p[2]) /3.0; 
c— t * t * t ; 
d— b*b— c ; 
i r (d>=0.0) 
{ 

d=pow ( ( sqrt (d) + fabs (b)) ,1.0/3.0) ; 

il (d! = 0.0) 

{ 
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3242 
3243 
3244 
3245 
3246 
3247 
3248 
3249 
3250 
3251 
3252 
3253 
3254 
3255 
3256 
3257 
3258 
3259 
3260 
3261 
3262 
3263 
3264 
3265 
3266 
3267 
3268 
3269 
3270 
3271 
3272 
3273 
3274 
3275 
3276 
3277 
3278 
3279 
3280 
3281 
3282 
3283 
3284 
3285 
3286 
3287 
3288 
3289 
3290 
3291 
3292 
3293 
3294 
3295 
3296 
3297 
3298 
3299 
3300 
3301 
3302 
3303 
3304 
3305 
3306 
3307 
3308 
3309 
3310 
3311 
3312 
3313 
3314 



if (b>0.0) 
b=-d; 

else 

b=d; 
c=t /b ; 



} 



d=r [2] [2] = sqrt\eqref {GrindEQ..0.75.}*(b-c) ; 
b=b+c ; 

c=r [1] [2]= -0.5*b-s ; 

if((b>0.0 \&\& s<=0.0) \textbar \toxtbar (b<0.0 \&\& s>O.Q)) 
{ 

r[l][l] = c; 
r[2][l]=-d; 
r[l][3] = b-s; 
r [ 2 ] [ 3 ] = . ; 

} 



r[l][l] = b-s; 
r [2] [1]=0.0; 
r[l][3] = c; 
r[2][3]=-d; 



} 

} /=t= end 2 equal or complex roots * / 
else 

{ 

if (b = = 0.0) 

d=atan \ eqrof { GrindEQ..1.0. } / 1 . 5 ; 

else 

d=atan ( sqrt (-d) /fabs (b) ) /3.0; 
if (b<0.0) 

b = 2.0* sqrt ( t ) ; 

else 

b=-2.0* sqrt ( t ) ; 
c— c o s ( d ) * b ; 

t=— sqrt \ eqref { Gr ind E Q __0_75 _ } * sin (d)*b — 0.5*c ; 
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3315 d=-t-c-s ; 

3316 

3317 c=c — s ; 

3318 

3319 t=t-s ; 

3320 

3321 i f ( f abs ( c ) > fabs ( t ) ) 

3322 

3323 { 
3324 

3325 r [1] [3] = c ; 

3326 

3327 } 
3328 

3329 else 
3330 

3331 { 
3332 

3333 r [ 1 ] [ 3 ] = t ; 

3334 

3335 t=c ; 

3336 

3337 } 
3338 

3339 i f ( fabs (d) > fabs ( t ) ) 

3340 

3341 { 
3342 

3343 r [ 1 ] [ 2 ] = d ; 

3344 

3345 } 
3346 

3347 else 
3348 

3349 { 
3350 

3351 r[l][2] = t; 

3352 

3353 t=d ; 

3354 

3355 } 
3356 

3357 r[l][l] = t; 

3358 

3359 for (k=l;k<4;k++) 

3360 

3361 r [2] [k]=0.0; 

3362 

3363 } 
3364 

3365 return ; 

3366 

3367 } 

3368 

3369 

3370 

3371 void BIQUADROOTS( double p [], double r[][5]) 

3372 

3373 { 

3374 

3375 /* 
3376 

3377 Array r[3][5] p [ 5 ] 

3378 

3379 Roots of poly p [0]* x\'{}4 + p [1]* x\'{}3 + p[2l*x\'{}2 + p[3]*x + p[4] = 



3380 

3381 x=r[l][k] + i r[2][k] k=l,...,4 

3382 

3383 * / 

3384 

3385 double a,b,c,d,o; 

3386 
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3387 


i n t k , j ; 




3388 






3389 


if(p[0] 


1.0) 


3390 






3391 


{ 




3392 






3393 


for (k- 


a;k<5;k++) 


3394 






3395 


p[k]^p[k]/p[0]; 


3396 






3397 


p[0]-l 


.0; 


3398 






3399 


} 




3400 






3401 


c^0.25*p [ 1 


]; 


3402 






3403 


b^2.0*e ; 




3404 






3405 


c— b* b ; 




3406 






3407 


c ; 




3408 






3409 


b-p[3] + b*( 


c-p[2]) ; 


3410 






3411 


a^p [2] - d ; 




3412 






3413 


c-p[4] + c^( 


c*a— p [3] ) ; 


3414 






3415 


a— a— d ; 




3416 






3417 


p[l]=0.5*a 




3418 






3419 


p[2] = (p[l]*p[l]-c) *0.25; 


3420 






3421 


p[3] = b*b/( 


-64.0) ; 


3422 






3423 


if (p[3]<0. 


0) 


3424 






3425 


{ 




3426 






3427 


CUBICROOTS ( p , r ) ; 


3428 






3429 


for (k = 


a;k<4;k++) 


3430 






3431 


{ 




3432 






3433 


if (r [2] [k]= = 0.0 \&\& r [i; 


3434 






3435 


{ 




3436 






3437 




d=r [ 1 ] [ k ] * 4 . ; 


3438 






3439 




a— a+d ; 


3440 






3441 




if (a>=0.0 \&\& b>=0. 


3442 






3443 




p[l]=sqrt(d) : 


3444 






3445 




else if (a< = 0.0 \&\& 


3446 






3447 




p[l]=sqrt(d) ; 


3448 






3449 




else 


3450 






3451 




p[l] = -sqrt (d) ; 


3452 






3453 




b = 0.5*(a+b/p [ 1] ) ; 


3454 






3455 




goto QUAD; 


3456 






3457 


} 




3458 






3459 


} 
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3460 

3461 } 
3462 

3463 i f ( p [ 2 ] < . ) 

3464 

3465 { 
3466 

3467 b=sqrt ( c ) ; 

3468 

3469 d=b+b-a ; 

3470 

3471 p[l]=0.0; 
3472 

3473 ir(d>0.0) 
3474 

3475 P[l]= sqrt (d) ; 

3476 

3477 } 
3478 

3479 else 
3480 

3481 { 
3482 

3483 if(p[l]>0.0) 
3484 

3485 b=sqrt (p [ 2] ) *2.0 + p [ 1] ; 

3486 

3487 c 1 s c 

3488 

3489 b=-sqrt (p [2] ) *2.0 + p [ 1 ] ; 

3490 

3491 if(b!=0.0) 
3492 

3493 { 
3494 

3495 p [ 1 ] = . ; 

3496 

3497 } 
3498 

3499 G 1 s o 

3500 

3501 { 
3502 

3503 for (k = l;k<5;k++) 

3504 

3505 { 
3506 

3507 r [ 1 ] [ k]=-c ; 

3508 

3509 r [2] [k]=0.0; 

3510 

3511 } 
3512 

3513 goto END; 

3514 

3515 } 
3516 

3517 } 
3518 

3519 QUAD: 
3520 

3521 p[2] = c/b; 

3522 

3523 QUADROOre ( p , r ) ; 

3524 

3525 lor (k=l;k<3;k++) 

3526 

3527 for(j=l;j<3;j++) 
3528 

3529 r [j ] [k+2] = r [j 1 [k] ; 

3530 

3531 p[i] = -p[i]; 

3532 
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3533 p[2] = b; 

3534 

3535 QUADROOTS ( p , r ) ; 

3536 

3537 loi (k=l;k<5;k++) 

3538 

3539 r [1] [k] = r [1] [k]-c; 

3540 

3541 END: 
3542 

3543 return; 

3544 

3545 } 



Listing 15. C-SOURCE CODE FOR UTILITY FUNCTIONS 



7. APPENDIX D. 

3547 program \ _c o n s t an t s . h : 

3548 
3549 
3550 

3551 // 



3552 

3553 //== INCLUDE ANSI C SYSTEM HEADER FILES 



3554 

3555 // 



3556 

3557 ^include <niath.h> // — for calls to trig , sqrt and power functions 

3558 

3559 

3560 

3561 // 



3562 

3563 //== DEFINE PROGRAM CONSTANTS 



3564 

3565 // 



3566 










3567 


#do 


f i n c 


NORMAL.TERMINATION 





3568 










3569 


#dc 


fine 


NOJNTERSECTION.POINTS 


100 


3570 










3571 


#dc 


fine 


ONEJNTERSECTIONJOINT 


101 


3572 










3573 


#dc 


fine 


LINE.TANGENT.TO.ELLIPSE 


102 


3574 










3575 


#dc 


f i ri c 


DISJOINT.ELLIPSES 


103 


3576 










3577 


#dc 


f i n c 


ELLIPSE2.0UTSIDETANGENT_ELLIPSE1 


104 


3578 










3579 


c f i n c 


ELLIPSE2JNSIDETANGENT.ELLIPSE1 


105 


3580 










3581 


#do 


fine 


ELLIPSES.INTERSECT 


106 


3582 










3583 


c f i n c 


TWOJNTERSECTION.POINTS 


107 


3584 










3585 


#dc 


fine 


THREEJNTERSECTION JOINTS 


108 


3586 










3587 


#do 


fine 


FOURJNTERSECTION JOINTS 


109 


3588 










3589 


#dc 


fine 


ELLIPSE1.INSIDE.ELLIPSE2 


110 


3590 










3591 


#dc 


f i n e 


ELLIPSE2.INSIDE.ELLIPSE1 


111 
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3592 












3593 


#dcl 


1 n e 


ELLIPSES_ARE_IDENTICAL 


112 




3594 












3595 


#del 


i 11 e 


INTERSECTION_POINT 


113 




3596 












3597 


#dcf 


i 11 e 


TANGENT_POINT 


114 




3598 












3599 












3600 












3601 


#doi 


i 11 e 


ERROR-ELLIPSEJARAMETERS 


-100 




3602 












3603 


#dcl 


i 11 e 


ERROR_DEGENERATE_ELLIPSE 


-101 




3604 












3605 


#dci 


i 11 e 


ERROR_POINTSJSIOT_ON_ELLIPSE 


— 102 




3606 












3607 


#dci 


i 11 e 


ERRORJNVERSE.TRIG 


-103 




3608 












3609 


#dci 


i 11 e 


ERRORiINE_POINTS 


— 104 




3610 












3611 


#dci 


i 11 e 


ERROR-QUARTIC.CASE 


-105 




3612 












3613 


#dcl 


i n e 


ERROR-POLYNOIvDALJJEGREE 


-107 




3614 












3615 


#dcl 


i n e 


ERROR_POLYNOMIAL_ROOTS 


— 108 




3616 












3617 


#dcf 


i 11 e 


ERRORJNTERSECTION_PTS 


-109 




3618 












3619 


#dcl 


i 11 e 


ERROR.CALCULATIONS 


— 112 




3620 












3621 












3622 












3623 


#dol 


i n e 


EPS 


+ 1.0E-07 




3624 












3625 


#doi 


i 11 e 


pi (2.0* asin (1.0) ) // — 


a maximum^ precision 


value of 


3626 












3627 


# d c f i n c 


twopi (2.0*pi) // — 


a maximmrtr- precision 


V alue of 


3628 












3629 












3630 












3631 












3632 












3633 












3634 












3635 


call 


_e s 


. c : 






3636 












3637 












3638 












3639 


#in elude <stdio.h> 






3640 












3641 


#inelude <math.h> 






3642 












3643 


#in elude " program_constants . h" 






3644 












3645 


d o u 1 


le 


ellipsc_segincnt (double A, double B, double XI, double Yl , 



3646 
3647 
3648 
3649 
3650 
3651 
3652 
3653 
3654 
3655 
3656 
3657 
3658 
3659 
3660 
3661 
3662 



double Y2 , int * MessageCodc ) 



i n t main ( i ii i argc , char ** argv) 
f 

double A, B; 
double XI, Yl; 
double X2, Y2 ; 
double arcal , area2 ; 
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3663 double pi = 2.0 * asin e q r e f { Gr indEQ 1_0_}; // — a maximum^ p r e ci s i o n 

value of pi 

3664 

3665 iiil rtn; 

3666 

3667 char msg [ 1 2 4 ] ; 

3668 

3669 printf ("Calling cllipsc_scgmcnt.ctcxtbackslasli n"); 

3670 

3671 

3672 

3673 // — case shown in Fig. 1 
3674 

3675 A = 4 . ; 

3676 

3677 B = 2 . ; 

3678 

3679 XI = 4./sqrt (5.); 

3680 

3681 Yl = 4./sqrt (5.); 

3682 

3683 X2 = -3.; 

3684 

3685 Y2 = -sqrt ( 7 . ) / 2 . ; 

3686 

3687 

3688 

3689 aroal = c 1 1 i p s o _s c g m o n t (A, B, XI, Yl , X2 , Y2 , &rtn); 

3690 

3691 sprintf (msg, "Fig 1: segment area — %15.8f , return_value — %d\ 

textbaekslash n", areal, rtn); 

3692 

3693 printf (msg) ; 

3694 

3695 

3696 

3697 // — case shown in Fig. 2 
3698 

3699 A = 4 . ; 

3700 

3701 B = 2 . ; 

3702 

3703 XI = -3.; 

3704 

3705 Yl = -sqrt ( 7 . ) / 2 . ; 

3706 

3707 X2 = 4./ sqrt (5. ) ; 

3708 

3709 Y2 = 4. / sqrt ( 5 . ) ; 

3710 

3711 

3712 

3713 area2 = e 1 1 i p s e _s e g m e n t (A, B, XI, Yl , X2 , Y2 , &rtn); 

3714 

3715 sprintf (msg, "Fig 2: segment area — %15.8f , return_value — % 

dtextbaekslash n", area2, rtn); 

3716 

3717 printf (msg) ; 

3718 

3719 

3720 

3721 sprintf (msg, "sum of ellipse segments — % 15 . 8 f t e x t b a c ks 1 as h n" , areal + 

arca2 ) ; 

3722 

3723 printf (msg) ; 

3724 

3725 sprintf ( msg ," total ellipse area by pi*a*b — %15.8ftextbaekslash n", pi* 

A*B) ; 

3726 

3727 printf (msg) ; 

3728 

3729 

3730 
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return rtn : 



3731 
3732 
3733 
3734 
3735 
3736 
3737 
3738 
3739 
3740 
3741 
3742 
3743 
3744 
3745 
3746 
3747 
3748 
3749 

3750 
3751 
3752 
3753 
3754 

3755 double \ t ex t b f { c II i p s c _ 1 i n c _ o v c r I a p } (double PHI, double A, double B 
double H , 

3756 
3757 
3758 
3759 
3760 
3761 
3762 

3763 int \ t e X t b f { main } (int argc , char argv ) 

3764 
3765 { 
3766 



call_el .c: 

^include <stdio.h> 
#include <math.h> 
#in elude " program_eonstants . h" 

double \textbf{cllipsc_segment} (double A, double B, double XI, double Yl , 
double X2, 

double Y2 , int * MessageCode ) ; 



double K, double XI, double Yl , double X2 , 
double Y2 , inL * MessagcCode ) ; 



3767 


double 


A, B; 


3768 






3769 


double 


H, K, PHI; 


3770 






3771 


double 


XI, Yl; 


3772 






3773 


double 


X2, Y2; 


3774 






3775 


double 


areal , arca2 ; 


3776 






3777 


double 


pi — 2.0 * \ t ex t b f { as i n } 






precision value of pi 



3778 
3779 
3780 
3781 
3782 
3783 
3784 
3785 
3786 
3787 
3788 
3789 
3790 
3791 
3792 
3793 
3794 
3795 
3796 
3797 
3798 
3799 
3800 



i n 1 rtn; 

c h a r msg [10 2 4]; 

\textbf{printf} ("Calling cllipse_line_ovcrlap.e\tcxtbackslash n"); 

// — case shown in Fig. 
A = 4. ; 
B = 2 . ; 
H = -6; 
K = 3; 

PHI = 3.* pi / 8.0; 
XI = -3.; 
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3801 Yl = 3 . ; 
3802 

3803 X2 = -7.; 

3804 

3805 Y2 = 7 . ; 

3806 

3807 

3808 

3809 arcal = \ t c x t b f { c 1 1 i p s c \ . 1 i n c \ .o v c r 1 a p } (PHI, A, B, H, K, XI, Yl , X2 , Y2 

, \&rtn ) ; 

3810 

3811 \ text bf { s pr i n t f } (msg,"Fig 4: area = \%15.8f , return.valuo = \%d\ 

textbaekslash n" , a real , rtn) ; 

3812 

3813 \ textbf { pr int f } (msg); 

3814 

3815 

3816 

3817 // — case shown in Fig. , points reversed 
3818 

3819 A = 4 . ; 

3820 

3821 B = 2 . ; 

3822 

3823 H = -6; 

3824 

3825 K = 3; 

3826 

3827 PHI = 3.* pi/ 8.0; 

3828 

3829 XI = -7.; 

3830 

3831 Yl = 7 . ; 

3832 

3833 X2 = -3.: 

3834 

3835 Y2 = 3 . ; 

3836 

3837 

3838 

3839 aroa2 = \ t ex t b f { c 1 1 i p s e \ . 1 i n e \ .o v e r 1 a p } (PHI, A, B, H, K, XI, Yl , X2 , Y2 

, \&rtn) ; 

3840 

3841 \ t c x t b f { s p r i n t f } (msg, "Fig 4 reverse; area — %15.8f , return_value — \%d\ 

textbaekslash n" , area2 , rtn) ; 

3842 

3843 \textbf{printf} (msg): 

3844 

3845 

3846 

3847 \ t ex t b f { s p r i n t f } (msg, "sum of ellipse segments — % 15 . 8 f t e x t b ac k s 1 as h n" , 

areal + area2 ) ; 

3848 

3849 \textbf{printf} (msg): 

3850 

3851 \ textbf { s pr i nt f} (msg 

ftextbaekslash n 

3852 

3853 \textbf{printf} (msg) 

3854 
3855 
3856 

3857 return rtn ; 

3858 
3859 } 
3860 
3861 
3862 
3863 
3864 

3865 call.ee . e : 
3866 
3867 



'total ellipse area by pi*a*b — %15. 
, pi *A*B) ; 
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3868 
3869 
3870 
3871 
3872 
3873 
3874 
3875 
3876 
3877 

3878 
3879 
3880 
3881 
3882 
3883 
3884 
3885 
3886 
3887 
3888 
3889 
3890 
3891 
3892 
3893 
3894 
3895 
3896 
3897 

3898 
3899 
3900 
3901 
3902 
3903 
3904 
3905 
3906 
3907 
3908 
3909 
3910 
3911 

3912 
3913 
3914 
3915 

3916 
3917 
3918 
3919 
3920 
3921 
3922 
3923 
3924 
3925 
3926 
3927 
3928 
3929 
3930 
3931 

3932 
3933 
3934 



^^include <stdio.h> 

:^in elude " program _eonstants . h" 

double e 1 1 i p s e _ c 1 1 i p s e _ o V e r 1 a p (double PHI_1 , double Al , double Bl , 

double HI, double Kl , double PHI.2 , 
double A2, double B2 , double H2 , double K2 

iut *rtnCodc) ; 



iut main ( j ]i t argc , eh 



a r * * a r 



{ 



double Al, Bl, HI, Kl , PHI_1 ; 

double A2, B2 , H2 , K2 , PHI_2 ; 

double area; 

i 11 1 r t n ; 

char msg [10 2 4] ; 

printf ("Calling cllipsc_ellipsc_ovcrlap.c\tcxtbackslash n\ 
t e xt b ac ks 1 ash n" ) ; 

/ / case 0—1 

Al ^ 3.; Bl ^ 2.; HI ^ 0.; Kl ^ 0.; PHI_1 ^ 0.; 

A2 ^ 2.; B2 ^ 1.; H2 ^ -.75; K2 ^ 0.25; PHI_2 ^ pi/4.; 

area — cllipsc_cllipsc_ovcrlap (PHI_1, Al, Bl, HI, Kl, 

PHI_2 , A2, B2, H2, K2 , \&;rtn); 

sprintf ( msg ,"Casc — 1: area — \%15.8f, return _value — \%d\ 
tcxtbackslash n'' , area , rtn) ; 

printf ( msg ) ; 

sprintf ( msg ," ellipse 2 area by pi*a2*b2 — \%15.Sf\ 

tcxtbackslash n" , pi*A2*B2) ; 

printf ( msg ) ; 
// case 0-2 

Al ^ 2.; Bl ^ 1.; HI ^ 0.; Kl ^ 0.: PHI_1 ^ 0.; 

A2 ^ 3.; B2 ^ 2.; H2 ^ -.3; K2 ^ -.25; PHI_2 ^ pi/4.; 

area — cllipse_cllipsc_ovcrlap (PHI_1, Al, Bl, HI, Kl, 

PHI_2 , A2, B2, H2, K2 , &:rtn}; 

sprintf ( msg ,"Case — 2: area — %15.8f, return\_value — \%d\ 
tcxtbackslash n'' , area , rtn) ; 

printf ( msg ) ; 
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3935 


sprinti (msg, ellipse 1 area by pi*al*bl = \%15.8i\ 




t ext b ac kslas h n" , pi*Al*Bl) ; 


3936 




3937 


• 

p r 1 n t f ( msg ) ; 


3938 




3939 




3940 




3941 


// case U—o 


3942 




3943 


Al — 2 . ; Bl — 1 . ; HI — . ; Kl — . ; PHI_1 — . ; 


3944 




3945 


An 1 1^. TUT r\ T . u o or;. T*" o i . td tu to -.-^ ; / a . 
Az — l.o; rSz — U.fo; ilz — — z , o ; rvz — i.o; r^rli_z — pi/4.; 


3946 




3947 


A TT 

area — ellipse_ellipse_overlap (PHI_1, Al , Bl , HI , Kl , 


3948 




3949 


DUT O AO "DO TUO T^O f7-T.4-T. ^ . 

r^rll_z , AZ , dZ , rlz , J\.z , &5rt n ) , 


3950 




3951 


sprinti (msg. Case — 3: area = \%15.8f , returii_value = \%d\ 




t cx t b ac k s 1 as h n" , area , rtn) ; 


3952 




3953 


■ 

p r 1 n 1 1 ( msg ) ; 


3954 




3955 


prmti ( lljllipses arc disjoint, ovelap area — U.U\ 




textbackslash n\textbackslash n" ) ; 


3956 




3957 




3958 




3959 


/ / c a s e 1—1 


3960 




3961 


Al = o.; Bl — z.; rll = u.; Kl — U.; Plll\_l — u.; 


3962 




3963 


AO • "Ro 1 • Tio 1 no/if^onQOfinnoo- t^o n o • pwto t-,^ 

AZ — ^ ■ ) DZ — J- ■ ) xlZ — — l.UZftOZUyZDUUZZj IvZ — U.ZO, j:^rli_Z — pi 


3964 




3965 


area — eiiipse_eiiipBe_overiap rai_ j. j , ^ ^ > j iv j. j 


3966 




3967 


T>TTT O AO "DO TUrO T^O \ P J \ . 

r'm_2 , Az , d2 , ri2 , K.Z , \&rt n J ; 


3968 




3969 


sprinti ( msg , Case 1 — 1: area — \%15.8r, return\_value = \%a\ 




textbackslash n" , area , rtn) ; 


3970 




3971 


• 

p r 1 n t f ( msg ) ; 


3972 




3973 


sprinti (msg, ellipse 2 area by pi*a2*b2 = \%15.8i\ 




textbackslash n" , pi*A2*B2) ; 


3974 




3975 


• 

p r 1 n 1 1 ( msg ) ; 


3976 




3977 




3978 




3979 


// — case. 1—2 


3980 




3981 


Al = z . ; Bl = 1 . ; Ml = U . ; K.1 = U . ; Fral.l = U . ; 


3982 




3983 


AO OK. TDO 1 ft. VO OO. irO 1 . TD TJ T O /A 

AZ — O.O, dZ — l.o, rlZ — . Z Z , ±\.Z — -J-i Jrrli_Z — pi/4.. 


3984 




3985 


area = ellipse_ellipse_overlap ( PHI_1 , Al , Bl , HI , Kl , 


3986 




3987 


T> T T TO AO "DO T T O Ty O \P J \. 

Pnl_2 , Az , Bz , ri2 , 1S.Z , \&rt n ) ; 


3988 




3989 


sprinti ( msg , Case 1 — 2: area = \7ol5.8r, retum_value = \%a\ 




textbackslash n" , area , rtn) ; 


3990 




3991 


print f (msg) ; 


3992 




3993 


sprintf (msg," ellipse 1 area by pi*albl = \%15.8f\ 




textbackslash n" , pi*Al*Bl) ; 


3994 




3995 


print f (msg) ; 


3996 




3997 




3998 




3999 


// — case 1—3 


4000 
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4001 


Al — 2.: ril — 1.: Ill — u.; Kl — U.: Flll_l — u.; 


4002 




4003 


Az — l.o; dA — U.Yo; rlz — — z . U 1 79 bo 9 oUo 5 ; Kz — 1.25; Frll_z — pi 


4004 




4005 


11' 11' 1 / T T T 1 A 1 1~^ 1 T T 1 1 

area — ellipse_eliipse_overlap (Flll_l, Al , Bl , HI , Kl , 


4006 




4007 


PHI_2 , A2, B2 , H2 , K2 , \&rtn); 


4008 




4009 


sprintt (msg, Case 1—3: area = %15.8r, re t urn \ _v al ue = \%a\ 
t ex t b ac k slas h n" , area , rtn) ; 


4010 




4011 


■ 

p r 1 n 1 1 ( msg ) ; 


4012 




4013 


printf (" Ellipses are disjoint, ovelap area — 0.0\ 
text backslash n \ text backslash n); 


4014 




4015 




4016 




4017 


//— case 1 


4018 




4019 


A 1 Q ■ TU 1 O . TU 1 n . Tj" 1 f 1 . TD TT T 1 n . 

Al — o . ; rSi — Z . . HI — U . ; Ki — U . ; r^rll_l — U . ; 


4020 




4021 


A2 — z , z o ; Ijz — i . o ; xlz — U . ; Kz — U.; r^rll_z — pi/4.; 


4022 




4023 


area = ellipse_ellipse_overlap y Fhll_l , Al , Bl , Ml , Kl , 


4024 




4025 


TITTT i\ AO no TTO T^O \ P J_ ^ 

FH1_2 , A2 , B2 , H2 , Kz , \&rt n J ; 


4026 




4027 


sprintt ( msg , Casc z — 1: area — \7ol5.or, return_value = \%a\ 
t c X t b a c k s 1 a s h n" , area , rtn) ; 


4028 




4029 


• 

printf ( msg ) ; 


4030 




4031 


s p r i n 1 1 ( msg ,' ellipse2area by pi*a2*b2 = \%15.8f\ 
t cx t b ac k slas h n , pi*A2*B2) ; 


4032 




4033 


printf ( msg ) ; 


4034 




4035 




4036 




4037 


// — case 2—2 


4038 




4039 


A 1 O . TI> 1 1 . XT 1 r\ iy 1 t\ . Ti XT T 1 r\ . 

Al = z . ; Bl = 1 . ; Ml = U . ; Kl = U . ; r'Ml_l = U . ; 


4040 




4041 


A o o . T30 1 T . xjo n . x/"o n . d tj t o ^ i I a . 

Az = o . ; Bz = 1 . ( ; Mz = U . ; Kz = U . ; Jr rll_z = p 1 / 4 . ; 


4042 




4043 


„11; „ „ii; i„ (ft XT T 1 A 1 TU 1 XT 1 iy i 

area = ellipse_ellipse_overlap ( f Ml_l , Al , Bl , Ml , Kl , 


4044 




4045 


T~i T T T O AO T^O TTO T^O \ V J_ \ 

FM1_2 , Az , B2 , M2 , Kz , \&rt n ) ; 


4046 




4047 


sprintt (msg, Case 2 — 2; area = \%15.8f, return_value = \%a\ 
t ext b ac kslas h n" , area , rtn) ; 


4048 




4049 


• 

printf ( msg ) ; 


4050 




4051 


sprinti ( msg , cllipselarea by pi*albl = \%l5.8f\ 

J. il 11 1^ "Al, TU 1 \ 

t c X t b ac k s 1 as n n , pi*Al*r>l); 


4052 




4053 


printf (msg) ; 


4054 




4055 




4056 




4057 


// — case 2—3 


4058 




4059 


Al = 3.; Bl = 2.; HI = 0.; Kl = 0.; PHI.l =0.; 


4060 




4061 


A2 = 2.; B2 = 1.; H2 = -2.; K2 = -1.; PHI.2 = pi /4.; 


4062 




4063 


area = e 11 i p s e .e 11 i p s e .o ve r 1 ap (PHI.l, Al , Bl , HI, Kl , 


4064 




4065 


PHI.2, A2, B2, H2, K2 , \&rtn); 


4066 
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sprintf (insg,"Case 2 — 3: area — \%15.8f, retum\_valuc - 
t c X t b a c k s 1 a s h n\ t c xt b ac ks 1 as li n" , area , rtn) ; 


- \/oa\ 


4068 






4069 


p r i 11 1 f ( insg ) ; 




4070 






4071 






4072 






4073 


/ /— case o — 1 




4074 






4075 


rtl — "J-? -Dl — ^-7 -til — U., Iv 1 — U., l^lli_l — u., 




4076 






4077 


i\Z — i->^ — J--) H-i — -L'l Ivz — U.oO, r n. — pi/4i., 




4078 






4079 


aiea — Giiipse_ciiipsG_oveiiap (^i^iii_i, ai, IjI, -tiJ-j J^J-i 




4080 






4081 


T3 TT T 9 AO UO TV'O 


\6t r t n j 


4082 






4083 


sprinti \ msg , Oasc o — li area — \/ol5.oi, return\_valuc - 
textbackslash. ii'' , area , rtn) ; 


- \/oa\ 


4084 






4085 


p r i 11 1 f ( msg ) ; 




4086 






4087 






4088 






4089 


/ /— — case . i— ic 




4090 






4091 


rti — ■^■5 1j1 — J--) HI — '-'■I ivi — i^iii_i — "J., 




4092 






4093 


Pi. J. — z . z ; I3Z — i . ; rlZ — U . o ; ivZ — U.; r^rll_Z — pi/4.; 




4094 






4095 


area — eiiipse_eiiipsc_oveiiap (^i^iii_i , IjI, -tiJ-i ^j-j 




4096 






4097 


TZ> TT T AO R9 TTO T^'O 
r^rli_Z , /\z , rSZ , rlZ , iVz , 


\&c r t n j 


4098 






4099 


sprinti (msg, Oase 3 — 2: area — \70i5.81 , r e t ur n \ _v al u e - 
textbackslash ii\textbackslasli 11" , area , rtn) ; 


- \/oa\ 


4100 






4101 


p r i n t f ( msg ) \ 




4102 






4103 






4104 






4105 


/ / case 4 — 1 




4106 






4107 


Al Q. T31 9> TT1 n - I<^1 n. T3TTT1 n. 

rtl — "J ■ ? T3l — ^ ■ 7 ^1 — U . , J:S.i — . , r^lli_l — U ■ i 




4108 






4109 


A2 — 3.; B2 — 1.; H2 — 1.; K2 — —0.5; PHI_2 — pi/4.; 




4110 






4111 


area — cIlipse_elIipsc_overlap (PHI_1, Al, Bl, HI, Kl, 




4112 






4113 


PHI_2 , A2, 32, H2, K2 , 


\&rtn) 


4114 






4115 


sprintf ( msg ,"Case 4—1: area — \%15.8f, return_value — 
textbackslash n'' , area , rtn) ; 


\%d\ 


4116 






4117 


p r i n t f ( msg ) ; 




4118 






4119 






4120 






4121 


return rtn; 




4122 






4123 } 
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